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Abstract 

There  are  a  number  of  interesting  applications  where  modeling  elastic  and/or  vis¬ 
coelastic  materials  is  fundamental,  including  uses  in  civil  engineering,  the  food  industry, 
land  mine  detection  and  ultrasonic  imaging.  Here  we  provide  an  overview  of  the  sub¬ 
ject  for  both  elastic  and  viscoelastic  materials  in  order  to  understand  the  behavior  of 
these  materials.  We  begin  with  a  brief  introduction  of  some  basic  terminology  and 
relationships  in  continuum  mechanics,  and  a  review  of  equations  of  motion  in  a  con¬ 
tinuum  in  both  Lagrangian  and  Eulerian  forms.  To  complete  the  set  of  equations,  we 
then  proceed  to  present  and  discuss  a  number  of  specific  forms  for  the  constitutive 
relationships  between  stress  and  strain  proposed  in  the  literature  for  both  elastic  and 
viscoelastic  materials.  In  addition,  we  discuss  some  applications  for  these  constitutive 
equations.  Finally,  we  give  a  computational  example  describing  the  motion  of  soil  ex¬ 
periencing  dynamic  loading  by  incorporating  a  specific  form  of  constitutive  equation 
into  the  equation  of  motion. 
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ships. 
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1  Introduction 


Knowledge  of  the  field  of  continuum  mechanics  is  crucial  when  attempting  to  understand 
and  describe  the  behavior  of  materials  that  completely  fill  the  occupied  space  and  thus  act 
like  a  continuous  medium.  There  are  a  number  of  interesting  applications  where  modeling  of 
elastic  and  viscoelastic  materials  is  fundamental.  One  interest  in  particular  is  in  describing 
the  response  of  soil  which  experiences  some  sort  of  impact.  This  may  result  from  buildings 
falling  or  being  imploded,  or  even  an  intentionally  introduced  impact  as  part  of  land  mine 
detection  efforts  (see  [48,  50]).  The  chief  interest  is  in  determining  what  would  happen  to 
buried  objects  given  a  particular  surface  impact.  In  the  case  of  a  building  implosion,  there  are 
concerns  for  buried  infrastructure  such  as  tunnels,  pipes,  or  nearby  building  infrastructure. 
Investigations  can  be  carried  out  to  determine  the  likely  forces  on  these  buried  objects  to 
ensure  that  the  force  delivered  into  the  soil  will  not  damage  other  infrastructure.  When 
detecting  land  mines,  the  methodology  developed  in  the  papers  cited  uses  an  impact  on  the 
ground  to  create  Rayleigh  surface  waves  that  are  subsequently  changed  upon  interacting 
with  a  buried  mine;  this  change  in  wave  form  might  be  detected  through  electromagnetic 
or  acoustic  means.  Creating  a  model  that  accurately  describes  these  Rayleigh  waves  is  key 
to  modeling  and  understanding  the  buried  land  mine  situation.  In  both  of  these  examples, 
one  must  study  the  soil  properties,  determine  a  valid  constitutive  relationship  for  the  soil, 
and  verify  the  accuracy  of  the  model.  One  can  then  use  the  models  to  predict  the  results 
from  different  forces,  soil  properties,  etc.  Another  application  is  the  non-invasive  detection 
of  arterial  stenosis  (e.g.,  see  [1,  3,  11,  34,  46]).  In  this  study,  blockages  in  the  artery  create 
turbulence  in  the  blood  flow,  which  then  generates  an  acoustic  wave  with  a  normal  and  shear 
component.  The  acoustic  wave  propagates  through  the  chest  cavity  until  it  reaches  the  chest 
wall,  where  a  series  of  sensors  detect  the  acceleration  of  the  components  of  the  wave.  The 
data  from  the  sensors  can  then  be  used  to  quickly  determine  the  existence  and  perhaps  the 
location  of  the  blockages  in  the  artery.  This  technique  is  inexpensive  and  non-invasive.  For 
such  a  technology  to  be  feasible,  a  mathematical  model  that  describes  the  propagation  of 
the  acoustic  wave  from  the  stenosis  to  the  chest  wall  will  be  necessary  to  correctly  detect 
the  location  of  a  blockage. 

The  goal  of  this  paper  is  to  provide  a  brief  introduction  of  both  elastic  and  viscoelastic 
materials  for  those  researchers  with  little  or  no  previous  knowledge  on  continuum  mechanics 
but  who  are  interested  in  studying  the  mechanics  of  materials.  Since  we  are  concerned  with 
purely  mechanical  theories,  thermodynamic  variables  such  as  temperature  and  entropy  are 
excluded  in  this  paper.  For  more  information  on  this  aspect,  the  interested  reader  may  refer 
to  [54] .  In  addition,  the  materials  that  we  are  considering  are  simple  (for  example  the  stress 
at  a  given  material  point  depends  only  on  the  history  of  the  first  order  spatial  gradient  of  the 
deformation  in  a  small  neighborhood  of  the  material  point  and  not  on  higher  order  spatial 
gradients)  and  non-aging  (the  microscopic  changes  during  an  experiment  can  be  neglected 
in  the  basic  model).  Our  presentation  is  part  tutorial,  part  review  but  not  a  comprehensive 
survey  of  a  truly  enormous  research  literature.  We  rely  on  parts  of  the  standard  literature 
and  discuss  our  view  of  generally  accepted  concepts.  We  present  a  discussion  of  topics  we 
have  found  useful  over  the  past  several  decades;  hence  approximately  20%  of  references  are 
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work  from  our  group.  We  have  not  meant  to  ignore  major  applications  in  the  many  fine 
contributions  of  others;  rather  our  presentation  reflects  a  certain  level  of  comfort  in  writing 
about  efforts  on  which  we  have  detailed  knowledge  and  experience. 

The  paper  is  outlined  as  follows:  In  Section  2  some  basic  terminology  (such  as  stress  and 
strain)  and  relationships  (e.g.,  the  relationship  between  strain  and  displacement)  of  contin¬ 
uum  mechanics  are  briefly  described.  In  addition,  we  give  a  brief  overview  of  equations  that 
describe  the  general  motion  of  a  continuous  material  in  both  Enlerian  and  Lagrangian  forms. 
These  equations  describe  motion  utilizing  a  general  relationship  between  stress  and  strain  for 
materials.  We  then  must  introduce  constitutive  relationships  (Section  3)  in  order  to  quantify, 
for  a  particular  material,  the  relationship  between  stress  and  strain.  This  is  where  we  see 
elasticity  and  viscoelasticity  arise  in  response  to  the  fact  that  the  stress-strain  relationships 
can  be  linear  or  nonlinear,  and  because  various  materials  exhibit  differing  properties  under 
stresses  and  strains.  Since  the  constitutive  relationships  are  material-specific,  a  wide  array 
of  approaches  to  model  these  relationships  has  been  developed;  this  is  where  we  focus  much 
of  our  attention  in  this  paper.  Once  this  material-dependant  determination  is  complete,  we 
are  able  to  solve  the  resulting  set  of  equations  with  proper  boundary  and  initial  conditions 
and  can  then  describe  the  motion  of  the  medium  under  different  inputs  such  as  an  impact 
or  force  applied  to  the  material.  We  conclude  the  paper  in  Section  4  by  giving  an  example 
describing  the  motion  of  soil  experiencing  dynamic  loading.  The  approach  illustrated  can 
be  of  use  in  a  number  of  applied  problems;  for  example,  in  predictive  analyses  and  possible 
interrogation  attempts  in  location  of  buried  targets  or  in  assessing  damage  to  underground 
infrastructures  due  to  perturbing  shocks  (earthquakes,  surface  explosions,  etc.). 


2  Preliminary  Notions  and  Balance  Laws 


Throughout  this  review,  bold  letters  are  used  to  denote  vectors  unless  otherwise  indicated, 

|  •  |  is  used  to  denote  the  determinant  of  a  matrix,  and  5^  denotes  the  Kronecker  delta,  that  is, 
Sij  =  1  for  i  —  j  and  zero  otherwise.  For  convenience  of  presentation,  we  may  occasionally  use 
the  Einstein  notational  convention  (or  index  convention),  where  the  convention  is  as  follows: 
the  repetition  of  an  index  in  a  term  will  denote  a  summation  with  respect  to  that  index  over 
its  range.  For  example,  the  Einstein  representations  for  EE  CijkiEki  find.  EE  dijdxidxj 

k  l  i  j 

are  given  by  CijkiSki  and  Sijdxidxj ,  respectively.  In  the  Cartesian  coordinate  system  we  may 
denote  coordinate  axes  by  x,  y,  and  z  or  by  x±,  X2  and  x^,  depending  again  on  the  ease  of 
presentation.  Accordingly,  the  components  of  a  tensor  a  are  denoted  by  axx,  axy,  crxz,  etc., 
in  reference  to  coordinates  x,  y  and  z,  and  are  denoted  by  ,  i,j  =  1,2,3  in  reference  to 
coordinates  x\,  X2  and  x%. 

In  this  section  we  first  introduce  some  basic  terminology  and  relationships  used  in  continuum 
mechanics  and  then  give  a  review  on  some  fundamental  physical  laws  such  as  the  conservation 
of  mass,  and  equation  of  motions  in  both  Lagrangian  and  Eulerian  forms.  The  content 
of  this  section  is  a  summary  of  material  from  several  popular  mechanics  books  including 
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[23,  25,  35,  40], 


2.1  Preliminaries 

2.1.1  Kinematics:  Deformation  and  Motion 

An  external  force  applied  to  an  object  results  in  a  displacement,  and  the  displacement  of  a 
body  generally  has  two  components: 

•  A  rigid-body  displacement:  in  this  case  the  relative  displacement  between  particles  is 
zero,  i.e.,  the  shape  and  size  of  the  body  does  not  change. 

•  A  deformation:  here  there  is  a  relative  displacement  between  particles,  i.e,  the  shape 
and/or  the  size  are  changed.  There  are  two  formulations  describing  deformation: 

—  Finite  strain  theory  which  deals  with  deformations  in  which  both  rotations  and 
strains  are  arbitrarily  large.  For  example,  elastomers,  plastically-deforming  ma¬ 
terials  and  other  fluids  and  biological  soft  tissue  often  undergo  large  deformations 
and  may  also  require  viscoelastic  ideas  and  formulations. 

—  Infinitesimal  strain  theory  which  treats  infinitesimal  deformations  of  a  continuum 
body.  Many  materials  found  in  mechanical  and  civil  engineering  applications, 
such  as  concrete  and  steel,  typically  undergo  small  deformations. 

For  this  presentation,  we  shall  focus  on  deformations.  When  analyzing  the  deformation  or 
motion  of  solids,  or  the  flow  of  fluids,  it  is  traditional  (and  helpful)  to  describe  the  sequence 
or  evolution  of  configurations  throughout  time.  One  description  for  motion  is  made  in  terms 
of  the  material  or  fixed  referential  coordinates,  and  is  called  a  material  description  or  the 
Lagrangian  description.  The  other  description  for  motion  is  made  in  terms  of  the  spatial 
or  current  coordinates,  called  a  spatial  description  or  Eulerian  description.  An  intuitive 
comparison  of  these  two  descriptions  would  be  that  in  the  Eulerian  description  one  places 
the  coordinate  or  reference  system  for  motion  of  an  object  on  the  object  as  it  moves  through 
a  moving  fluid  (e.g.,  on  a  boat  in  a  river)  while  in  the  Lagrangian  description  one  observes 
and  describes  the  motion  of  the  object  from  a  fixed  vantage  point  (e.g.,  motion  of  the  boat 
from  a  fixed  point  on  a  bridge  over  the  river  or  on  the  side  of  the  river.). 


Lagrangian  Description  In  a  Lagrangian  description  an  observer  standing  in  the  refer¬ 
ential  frame  observes  the  changes  in  the  position  and  physical  properties  as  the  material 
particles  move  in  space  as  time  progresses.  In  other  words,  this  formulation  focuses  on  indi¬ 
vidual  particles  as  they  move  through  space  and  time.  This  description  is  normally  used  in 
solid  mechanics.  In  the  Lagrangian  description,  the  motion  of  a  continuum  is  expressed  by 
the  mapping  function  h  given  by 

x  =  h(X,  t),  (2.1) 
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which  is  a  mapping  from  initial  (undeformed/material)  configuration  Oo  to  the  present  (de¬ 
formed/spatial)  configuration  Qt-  Hence,  in  a  Lagrangian  coordinate  system  the  velocity  of 
a  particle  at  X  at  time  t  is  given  by 

<9x  <9h(X, 

)t)  =  ~dt=  dt 

and  the  total  derivative  (or  material  derivative)  of  a  function  which  is  denoted  by  a 

dot  or  the  symbol  — ,  is  just  the  partial  derivative  of  with  respect  to  t, 


D 

Dt 


V>(X,t) 


d_ 

dt 


^(X,t). 


Eulerian  Description  An  Eulerian  description  focuses  on  the  current  configuration  Dt, 
giving  attention  to  what  is  occurring  at  a  moving  material  point  in  space  as  time  progresses. 
The  coordinate  system  is  relative  to  a  moving  point  in  the  body  and  hence  is  a  moving 
coordinate  system.  This  approach  is  often  applied  in  the  study  of  fluid  mechanics.  Math¬ 
ematically,  the  motion  of  a  continuum  using  the  Eulerian  description  is  expressed  by  the 
mapping  function 

X  =  h-1(x,  t), 

which  provides  a  tracing  of  the  particle  which  now  occupies  the  position  x  in  the  current 
configuration  from  its  original  position  X  in  the  initial  configuration  Ho-  The  velocity  of 
a  particle  at  x  at  time  t  in  Eulerian  coordinate  system  is 

v(x,t)  =  V(h~1(x,f),t). 

Hence,  in  an  Eulerian  coordinate  system  the  total  derivative  (or  material  derivative)  of  a 
function  if(x,t)  is  given  by 

D  d  3 '  Q  Q 

— r/;(x,  t)  =  t^(x,  t)  +  *)  =  dt ^X’  ^  +  V(X’  ^  '  V^(x’ 

i=  1  1 

Remark  2.1.  There  are  a  number  of  the  different  names  often  used  in  the  literature  to 
refer  to  Lagrangian  and  Eulerian  configurations.  Synonymous  terminology  includes  ini¬ 
tial/referential,  material,  undeformed,  fixed  coordinates  for  Lagrangian  and  current /present, 
space,  deformed,  moving  coordinates  for  Eulerian  reference  frames. 


2.1.2  Displacement  and  Strain 

A  particle  P  located  originally  at  the  coordinate  X  =  (X1:  X2,  X3)T  is  moved  to  a  place  P' 
with  coordinate  x  =  (x±,  x2,  X3)T  when  the  body  moves  and  deforms.  Then  the  vector  PP7. 
is  called  the  displacement  or  deformation  vector  of  the  particle.  The  displacement  vector  is 

x-X.  (2.2) 


5 


Let  the  variable  X  =  (Xi,X2,X3)t  identify  a  particle  in  the  original  configuration  of  the 
body,  and  x  =  (xi,  X2,  x3)T  be  the  coordinates  of  that  particle  when  the  body  is  deformed. 
Then  the  deformation  of  a  body  is  known  if  xll  X2  and  x3  are  known  functions  of  Ah,  X2, 

X3: 

Xi  =  Xi(X  1,X2,X3),  i  =  1,2,3. 

The  (Lagrangian)  displacement  of  the  particle  relative  to  X  is  given  by 

U(X)  =  x(X)  -  X.  (2.3) 

If  we  assume  the  transformation  has  a  unique  inverse,  then  we  have 

Xi  =  Xi{x  1,  x2,  x3),  i  =  1,  2,  3, 

for  every  particle  in  the  body.  Thus  the  (Eulerian)  displacement  of  the  particle  relative  to 
x  is  given  by 

u(x)  =  x  —  X(x).  (2.4) 

To  relate  deformation  with  stress,  we  must  consider  the  stretching  and  distortion  of  the  body. 
For  this  purpose,  it  is  sufficient  if  we  know  the  change  of  distance  between  any  arbitrary  pair 
of  points. 

Consider  an  infinitesimal  line  segment  connecting  the  point  P( Xi,X2,X3)  to  a  neighboring 
point  Q(X  1  +  dXi,X2  +  dX2,X3  +  cLA3)  (see  Figure  1).  The  square  of  the  length  of  PQ  in 


Figure  1:  Deformation  of  a  body. 

the  original  configuration  is  given  by 

|dX|2  =  (dX)TdX  =  (dX3)2  +  (dX2)2  +  (dX3f. 

When  P  and  Q  are  deformed  to  the  points  P\x  1,  X2,  x3)  and  Q'(x  1  +  dx i,x2  +  dx2,  x3  +  dx 3), 
respectively,  the  square  of  length  of  P'Q'  is 

|dx|2  =  (dx)Tdx  =  (dx  1)2  +  (dx2)2  +  (dx3)2 . 
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Definition  2.2.  The  configuration  gradient  (often,  in  something  of  a  misnomer,  referred  to 
as  deformation  gradient  in  the  literature)  is  defined  by 


dxi 

dx\ 

dxi 

dXx 

dX2 

dX3 

dx 

dx2 

dx2 

dx2 

dX  = 

dXx 

dX2 

dx3 

dx3 

dx3 

dx3 

L  0Xx 

dX2 

dx3  J 

The  Lagrangian  strain  tensor  The  Lagrangian  strain  tensor  is  measured  with  respect 
to  the  initial  configuration  (i.e. ,  Lagrangian  description).  By  the  definition  of  configuration 
gradient,  we  have  dx  =  AdX  and 

|dx|2  —  |dX|2  =  (dx)Tdx  —  (dX)TdX 

=  (dX)TATAdX  -  (dX)TdX 
=  {d,X)T(ATA  —  J)dX. 

The  Lagrarigian  ( finite )  strain  tensor  E  is  defined  by 

E  =  i(.4T,4-/).  (2.6) 

The  strain  tensor  E  was  introduced  by  Green  and  St.  Venant.  Accordingly,  in  the  literature 
E  is  often  called  the  Green’s  strain  tensor  or  the  Green- St.  Venant  strain  tensor.  In  addition, 
from  (2.6)  we  see  that  the  Lagrangian  strain  tensor  E  is  symmetric. 


If  the  strain  satisfies  ATA  —  1  =  0,  then  we  say  the  object  is  undeformed;  otherwise,  it  is 
deformed.  We  next  explore  the  relationship  between  the  displacement  and  strain.  By  (2.3) 
and  (2.5)  we  have  the  true  deformation  gradient  given  by 


VU  = 


or 


dxi 

~dXl  ~ 

dx2 


dx\ 


dX2 
dx2 


dU.j 


dXx 

dx3 

dXx 

dx; 


dx2 

dx3 


-  1 


OX, 


dxi 

dX~s 
dx2 
dX~3 
dx3  _ 
dX3 


=  A  —  I 


dX,  dX 


j  —  1)  2,  3,  i  —  1,  2,  3, 


'■j 

where  /  is  the  identity  matrix,  and  U  =  (fifi,  U2,  U3)T .  Thus,  because  A  =  VU  +  /,  the 
relationship  between  Lagrangian  strain  (2.6)  and  displacement  is  given  by 


E. 


v 


1  r  dUi  du3  duk  duk 

2  [dx j  +  dXt  +  dXi  dXj 


where  Eij  is  the  (i,j)  component  of  strain  tensor  E. 


7 


The  Eulerian  strain  tensor  The  Eulerian  strain  tensor  is  measured  with  respect  to  the 
deformed  or  current  configuration  (i.e.,  Eulerian  description).  By  using  dX  =  T-1dx,  we 
find 

|dx|2  —  |dX|2  =  (dx)Tdx  —  (dx)r(T_1)TT_1dx 
=  (dx)T(/  -  (^-1)T^-1)dx 

and  the  Eulerian  (finite)  strain  tensor  e  is  defined  by 


e  =  -jI-(A-')TA~l). 


(2.7) 


The  strain  tensor  e  was  introduced  by  Cauchy  for  infinitesimal  strains  and  by  Almansi  and 
Hamel  for  finite  strains;  e  is  also  known  as  Almansi’s  strain  in  the  literature.  In  addition, 
we  observe  from  (2.7)  that  Eulerian  strain  tensor  e  is  also  symmetric. 

We  also  may  give  the  relationship  between  the  displacement  and  strain  in  the  Eulerian 
formulation.  By  (2.4)  and  (2.5)  we  have 

_  dX] 
dx\ 

dX2 


Vu  = 


dxi 

dXs 

dxi 


dXl 

dXi  1 

dx2 

dx3 

dX2 

dX2 

dx2 

dx3 

dX3 

x  dX3 

dx2 

dx3  . 

-  I  -A 


-i 


or 


duj  .  dXj 

=  j  =  1,2,3,  i  =  1,2,3, 


where  u  =  (u\,  u2,  u3)T . 
given  by 


dxj  3  dxj 

Thus,  the  relationship  between  Eulerian  strain  and  displacement  is 


C-ij  —  — 
13  2 


dui  duj  duk  duk 
dxi  dx; 


dxi  dxj 


where  is  the  (i,j)  component  of  strain  tensor  e. 

Remark  2.3.  There  are  two  tensors  that  are  often  encountered  in  the  finite  strain  theory. 
One  is  right  Cauchy- Green  configuration  (deformation)  tensor ,  which  is  defined  by 


Dr  =  = 


dxk  dxk 
dX dX* 


and  the  other  is  the  left  Cauchy- Green  configuration  (deformation)  tensor  defined  by 

Dt  = AAT = 


dxi  dxj 
dXk  3Xk 


The  inverse  of  DL  is  called  the  Finger  deformation  tensor.  Invariants  of  Dr  and  DL  are 
often  used  in  the  expressions  for  strain  energy  density  functions  (to  be  discussed  below  in 
Section  3.1.2).  The  most  commonly  used  invariants  are  defined  to  be  the  coefficients  of  their 


characteristic  equations.  For  example,  frequently  encountered  invariants  of  DR  are  defined 
by 

h  —  tr(-D/?)  =  +  A2  +  Ag, 

^2  =  2  —  (M-Dr))2]  =  A1A2  +  A2A3  +  A3A3, 

h  =  det(DR)  =  XjXjXl 

where  A*,  i  =  1,2,3  are  the  eigenvalues  of  A,  and  also  known  as  principal  stretches  (these 
will  be  discussed  later  in  Section  3.1.2). 


Infinitesimal  strain  theory  Infinitesimal  strain  theory  is  also  called  small  deformation 
theory ,  small  displacement  theory  or  small  displacement- gradient  theory.  In  infinitesimal 
strain  theory,  it  is  assumed  that  the  components  of  displacement  tq  are  such  that  their  first 
derivatives  are  so  small  that  higher  order  terms  such  as  the  squares  and  the  products  of  the 
partial  derivatives  of  ul  are  negligible  compared  with  the  first-order  terms.  In  this  case  e^ 
reduces  to  Cauchy’s  infinitesimal  strain  tensor 


£ij  2 


dm  du,j 


+ 


dxn  dx 


(2.8) 


(hence  the  infinitesimal  strain  tensor  is  also  symmetric).  Thus,  =  e^.  We  note  that 
in  this  case  the  distinction  between  the  Lagrangian  and  Eulerian  strain  tensors  disappears 
(i.e.,  Eij  Ri  eij  ~  £^-),  since  it  is  immaterial  whether  the  derivatives  of  the  displacement  are 
calculated  at  the  position  of  a  point  before  or  after  deformation.  Hence,  the  necessity  of 
specifying  whether  the  strains  are  measured  with  respect  to  the  initial  configuration  (La¬ 
grangian  description)  or  with  respect  to  the  deformed  configuration  (Eulerian  description)  is 
characteristic  of  a  finite  strain  analysis  and  the  two  different  formulations  are  typically  not 
encountered  in  the  infinitesimal  theory. 


2.1.3  Stress 


Stress  is  a  measure  of  the  average  amount  of  force  exerted  per  unit  area  (in  units  N/m 2  or 
Pa),  and  it  is  a  reaction  to  external  forces  on  a  surface  of  a  body.  Stress  was  introduced  into 
the  theory  of  elasticity  by  Cauchy  almost  two  hundred  years  ago. 


Definition  2.4.  The  stress  vector  (traction)  is  defined  by 

dF  ,  AF 
T(  =  =  hm  — — , 

dT  Ar^o  AT  ' 


where  the  superscript  (n)  is  introduced  to  denote  the  direction  of  the  normal  vector  n  of  the 
surface  T  and  F  is  the  force  on  the  surface. 


To  further  elaborate  on  this  concept,  consider  a  small  cube  in  the  body  as  depicted  in  Figure 
2  (left).  Let  the  surface  of  the  cube  normal  (perpendicular)  to  the  axis  z  be  donated  by  Vrz. 
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Figure  2:  Notations  of  stress  components. 


Let  the  stress  vector  that  acts  on  the  surface  VTZ  be  T(e3\  where  e3  =  (0,0,1).  Resolve 
T(e3^  into  three  components  in  the  direction  of  the  coordinate  axes  and  denote  them  by  azx, 
azy  and  azz.  Similarly  we  may  consider  surface  VT^,  and  VT,y  perpendicular  to  x  and  y, 
the  stress  vectors  acting  on  them,  and  their  components  in  the  x,  y  and  z  directions.  The 
components  axx ,  ayy  and  <jzz  are  called  normal  stresses,  and  axy,  crxz,  cryx,  ayz,  azx  and  azy 
are  called  shear  stresses.  A  stress  component  is  positive  if  it  acts  in  the  positive  direction 
of  the  coordinate  axes.  We  remark  that  the  notation  <7face  cprcction  consistently  used  in 
elasticity  theory. 


Definition  2.5.  The  Cauchy  stress  tensor  is  defined  by 


'J'(ei) 

@ XX 

GXy 

GXz 

(7  = 

rJ'(e2) 

= 

&yx 

ayy 

Gyz 

rp(e3) 

O zx 

Gzy 

G  ZZ 

where  ei  =  (1,  0,  0),  e2  =  (0, 1,  0)  and  e3  =  (0,  0, 1). 


(2.9) 


We  have  the  following  basic  formulation  due  to  Cauchy. 

Theorem  2.6.  [24,  p.69]  Let  T(n)  be  the  stress  vector  acting  on  dT  whose  outer  normal 
vector  is  n,  as  illustrated  in  Figure  3.  Cauchy’s  formula  expresses  T(n)  as  a  function  of 
the  stress  vectors  on  the  planes  perpendicular  to  the  coordinate  axes,  i.e.,  in  terms  of  the 
components  of  the  Cauchy  stress  tensor.  This  formula  asserts  that 

Tj.  ''  Gxxnx  T  (Tyxny  T  Gzxnz , 

Ty  i  G xyl ~lx  +  Gyyny  T  Gzynz,  (^.10) 

T[  i  GxzTlx  T  Gy  z  fly  +  Gzznz. 


Here  T<n>  =  ,  and  n  =  (■ nx,ny,nz)T . 

written  concisely  as 


rj'(n)  _  -T 


a  n, 


where  a  is  the  Cauchy  stress  tensor  defined  in  (2.9). 


Cauchy’s  formula  (2.10)  can  be 
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j(n) 


Figure  3:  Stress  vector  acting  on  a  plane  with  normal  n. 

Remark  2.7.  In  addition  to  the  Cauchy  stress  tensor,  there  are  other  stress  tensors  encoun¬ 
tered  in  practice,  such  as  the  first  Piola-Kirchhoff  stress  tensor  and  second  Piola-Kirchhoff 
stress  tensor.  The  differences  between  the  Cauchy  stress  tensor  and  the  Piola-Kirchhoff 
stress  tensors  as  well  as  relationships  between  the  tensors  can  be  illustrated  as  follows: 


•  Cauchy  stress  tensor:  relates  forces  in  the  present  (deformed/spatial)  configuration  to 
areas  in  the  present  configuration.  Hence,  sometimes  the  Cauchy  stress  is  also  called 
the  true  stress.  In  addition,  the  Cauchy  stress  tensor  is  symmetric,  which  is  implied 
from  the  fact  that  the  equilibrium  of  an  element  requires  that  the  resultant  moment 
vanish.  We  will  see  in  Section  2.2.4  that  the  Cauchy  stress  tensor  is  used  in  the  Eulcrian 
equation  of  motion.  Hence  the  Cauchy  stress  tensor  is  also  referred  to  as  the  Eulcrian 
stress  tensor. 

•  First  Piola-Kirchhoff  stress  tensor  (also  called  the  Lagrangian  stress  tensor  in  [23,  24]): 
relates  forces  in  the  present  configuration  with  areas  in  the  initial  configuration.  The 
relationship  between  the  Erst  Piola-Kirchhoff  stress  tensor  P  and  the  Cauchy  stress 
tensor  cr  is  given  by 

P  =  \A\(tA~t  .  (2.11) 

From  the  above  equation,  we  see  that  in  general  the  first  Piola-Kirchhoff  stress  tensor 
is  not  symmetric  (its  transpose  is  called  the  nominal  stress  tensor  or  engineering  stress 
tensor).  Hence,  the  first  Piola-Kirchhoff  stress  tensor  will  be  inconvenient  to  use  in  a 
stress-strain  law  in  which  the  strain  tensor  is  always  symmetric.  In  addition,  we  will 
see  in  Section  2.2.5  that  the  first  Piola-Kirchhoff  stress  tensor  is  used  in  the  Lagrangian 
equation  of  motion.  As  pointed  out  in  [24] ,  the  first  Piola-Kirchhoff  stress  tensor  is  the 
most  convenient  for  the  reduction  of  laboratory  experimental  data. 

•  Second  Piola-Kirchhoff  stress  tensor  (referred  to  as  Kir choff  stress  tensor  in  [23,  24], 
though  in  some  references  such  as  [40]  Kirchhoff  stress  tensor  refers  to  a  weighted 
Cauchy  stress  tensor  and  is  defined  by  |A|cr):  relates  forces  in  the  initial  configura- 
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tion  to  areas  in  the  initial  configuration.  The  relationship  between  the  second  Piola- 
Kirchhoff  stress  tensor  S  and  the  Cauchy  stress  tensor  cr  is  given  by 

S  =  ^A-'crA-7.  (2.12) 

From  the  above  formula  we  see  that  the  second  Piola-Kirchhoff  stress  tensor  is  sym¬ 
metric.  Hence,  the  second  Piola-Kirchhoff  stress  tensor  is  more  suitable  than  the  first 
Piola-Kirchhoff  stress  tensor  to  use  in  a  stress-strain  law.  In  addition,  by  (2.11)  and 
(2.12)  we  find  that 

S  =  H_1P.  (2.13) 

Note  that  for  infinitesimal  deformations,  the  Cauchy  stress  tensor,  the  first  Piola-Kirchhoff 
stress  tensor  and  the  second  Piola-Kirchhoff  tensor  are  identical.  Hence,  it  is  necessary 
only  in  finite  strain  theory  to  specify  whether  the  stresses  are  measured  with  respect  to  the 
initial  configuration  (Lagrangian  description)  or  with  respect  to  the  deformed  configuration 
(Eulerian  description). 


2.2  Equations  of  Motion  of  a  Continuum 

There  are  two  common  approaches  in  the  literature  to  derive  the  equation  of  motion  of 
a  continuum:  the  differential  equation  approach  (for  example,  in  [25])  and  the  integral 
approach  (for  example,  in  [25,  40]).  In  this  section,  we  will  use  an  integral  approach  to 
derive  the  equation  of  continuity  and  the  equation  of  motion  of  a  continuum,  first  in  the 
Eulerian  (or  moving)  coordinate  system  and  then  in  the  Lagrangian  coordinate  system.  In 
the  following,  we  will  sometimes  use  to  denote  Qt,  and  Y  to  denote  Yt  for  ease  in  the 
presentation;  either  will  refer  to  a  volume  element  that  is  time  dependent. 

Before  deriving  the  equations  of  motion  of  a  continuum,  we  first  discuss  forces.  There  are 
two  types  of  external  forces  acting  on  material  bodies  in  the  mechanics  of  continuum  media: 


•  Body  forces  ( N/m 3),  acting  on  elements  of  volume  of  body.  For  example,  gravitational 
forces  and  electromagnetic  forces  are  body  forces. 

•  Surface  forces  (N/m2),  or  stress,  acting  on  surface  elements.  For  example,  aerodynam¬ 
ics  pressure  acting  on  a  body,  stress  between  one  part  of  a  body  on  another,  etc.,  are 
surface  forces. 


Then  the  total  force  F  acting  upon  the  material  occupying  the  region  H  interior  to  a  closed 
surface  Y  is 

F=  <f  T{n)dY  +  I  fdn,  (2.14) 

Jr  Jn 

where  f  =  (fx,  fy,  fz)T  is  the  body  force,  and  T(n)  is  the  stress  vector  acting  on  dY  whose 
outer  normal  vector  is  n. 
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The  expression  (2.14)  is  a  universal  force  balance  statement  independent  of  any  particular 
coordinate  system  being  used.  Of  course,  with  either  the  Eulerian  or  Lagrangian  formulation, 
the  stresses  and  forces  must  be  expressed  in  terms  of  the  appropriate  coordinate  system. 


2.2.1  The  Material  Derivative  of  a  Volume  Integral 


To  carry  out  our  derivations,  we  need  a  calculus  for  interchanging  integration  and  differentia¬ 
tion  when  both  the  limits  of  the  integration  and  the  integrand  depend  on  the  differentiation 
variable.  Let  <£>(£)  be  a  volume  integral  of  a  continuously  differential  function  4>(x,y,  z,t) 
defined  over  a  spatial  domain  Qt  occupied  by  a  given  set  of  material  particles  at  time  t, i.e., 


$(f)  =  ///  (j)(x,y,  z,t)dxdydz. 


Then  the  rate  of  change  of  <£>(£)  with  respect  to  t  is  given  by  (suppressing  the  multiple 
integral  notation  here  and  below  when  it  is  clearly  understood  that  the  integral  is  a  volume 
or  surface  integral) 


dd> 

dt 


'Fit 


<90 

~dt 


dtt  +  I  (4>vxnx  +  <\>vyny  +  (pvznz)dr, 


i  rt 


(2.15) 


where  on  the  boundary  Tt  of  Qt,  v  =  v(t)  is  the  velocity  v(f) 
(2.15)  can  be  written  concisely  as 


dd> 

dt 


'Fit 


<90 

~dt 


dQ+  0v  •  ndT. 


'Ft 


f dx  dy  dz\T 
\dt’  dt’  dt ) 


Equation 


The  first  term  on  the  right  side  corresponds  to  rate  of  change  in  a  fixed  volume,  and  the 
second  term  corresponds  to  the  convective  transfer  through  the  surface.  By  Gauss’  theorem, 
equation  (2.15)  can  also  be  written  as 


d$  _  f  fdcf)  d(f)vx  d(f)Vy  d(pvz\ 

Tt  -  V ai +  ~aT  +  “ay  ~aTj 


(2.16) 


or  more  concisely  as 


dt 


W+v-(*v])da 


This  rate,  called  the  material  derivative  of  <f>,  is  defined  for  a  given  set  of  material  particles 
in  a  moving  volume.  We  note  that  when  Qt  =  G0  f°r  all  t  (i.e.,  the  boundary  T  is  not  moving 
so  that  v  =  0),  this  becomes  simply 


d_ 

dt 


4>(x,  y,  z ,  t)dxdydz 


'Fl0 


<90 

~dt 


(. x ,  y ,  z,  t)dxdydz. 
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2.2.2  The  Equation  of  Continuity 


We  next  derive  the  equation  of  continuity  for  an  arbitrary  mass  of  particles  that  may  be 
changing  in  time.  The  mass  contained  in  a  domain  Qt  at  time  t  is 


m(t)  —  /  p(x,y,  z,t)dxdydz. 


Conservation  of  mass  requires  that 


dm 

dt 


0  and  thus  we  have  from  (2.16) 


dm 

dt 


dp_  d  pvx  dpvy  d  pvz 
dt  dx  dy  dz 


dQ. 


Hence,  we  obtain 


dp_  dpvx  dpVy  dpvz 
dt  dx  dy  dz 


0. 


Since  the  above  equality  holds  for  an  arbitrary  domain  Qt,  we  obtain  the  pointwise  equation 
of  continuity 

dff  dpvx  dpvy  dpvz 

dt  dx  dy  dz 

which  can  be  written  concisely  as 


0, 


(2.17) 


|  +  V.(pv)=0. 


2.2.3  Reynolds  Transport  Theorem 


In  this  subsection,  we  will  use  the  material  derivative  (2.16)  as  well  as  the  equation  of 
continuity  (2.17)  to  derive  the  celebrated  Reynolds  transport  theorem.  By  (2.16)  we  find 
that 

d  f  /'  (d(pvz)  dpvzvx  dpvzvy  dpvzvz\ 

*  L  = L  {— + — + — + —) dn- 

Then  using  the  equation  of  continuity  (2.17),  we  find  that  the  integrand  of  the  right  side  of 
the  above  equation  is  equal  to 


dp 


dv. 


mv‘  +  p~m+  ^ 


dpvx  dpv.y  dpvz 

dx  dy  dz 


dv. 


dv. 


dv. 


+  P^  +  P^  +  P^-J7 


'dp  dpvx  dpv.y  dpvz\  ( dv,  dvz  dvz  dvz 

-v^fr+lrr  +  ~ai  +  ~dr)+,,VW  +  Vx-fo+v>-^  +  v‘-fc 


f  dvz 


dv. 


dv. 


dv. 


=  P{-W  +  V^  +  V^  +  Vz^ 
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Hence,  we  have 


d 


f  dv 


dvz 


dvz 


dv. 


—  /  pvzdtt  =  p  —  +  vx—  +  vy—  +  vz—  dQ. 


dt 


'Clt 


1 nt 


V  dt 


dx 


dy 


dz 


(2.18) 


Equation  (2.18)  is  the  Reynolds  transport  theorem,  which  is  usually  written  concisely  as 


■j:  [  pvzdU  = 


dt 


P- 


Dv , 


'nt 


Dt 


dD, 


where 


Dvz 

Dt 


is  the  total  derivative  of  vz,  and  is  given  by 


D  dvz  dvz  dvz  dvz 

-vM,  y,  z,  t)  =  w  +  v^ 


We  note  that  the  above  is  independent  of  any  coordinate  system  and  depends  only  on  the 
rules  of  calculus  and  the  assumptions  of  continuity  of  mass  in  a  time  dependent  volume  of 
particles. 


2.2.4  The  Eulerian  Equations  of  Motion  of  a  Continuum 

We  are  now  ready  to  use  the  above  rules  of  calculus  and  the  continuity  of  mass  as  embodied  in 
the  Reynolds  transport  theorem  to  derive  the  equations  of  motion  in  an  Eulerian  coordinate 
system.  Throughout  we  have  0  =  and  T  =  Tt  (we  will  suppress  the  subscripts)  and  we 
assume  the  coordinate  system  (x,  y,  z)  is  now  moving  (changing  with  the  volume  element) 

with  a  velocity  v  =  ( — .  — .  —  |  of  the  deformation  of  the  material.  The  resultant  force 
V  dt  dt  dt ) 

Fz  in  the  z- direct  ion  on  an  arbitrary  volume  is 


Fz  =  /  T^dV  +  /  fzJm. 


(2.19) 


By  Cauchy’s  formula  (2.10)  and  Gauss’  theorem  we  have 

[  ^  dV  f  ( O'xzV'x  T  Gyz'R'y  T  (OzzTlz)  dT 


daxz  davz  dazz  . 

-r^-  +  ~~zr—  +  1  dfi. 

in  \  dx  dy  dz 

Hence,  by  the  above  equality  and  (2.19)  we  obtain 


Fz  = 


dcrxz  d(jyZ  duzz  . 

+  fz  dn. 

ox  ay  oz 


Newton’s  law  states  that 

d  i 


pMf!=  t  (f^  +  ^  +  f^  +  Z)dn. 

dt  Jfi  Jn  \  dx  dy  dz 
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Hence,  by  the  Reynolds  transport  theorem  we  have 


dvz 

dt 


+  vx 


dvz 

dx 


+  v, 


dv? 


dv? 


Z  ^UZ 

y  dyVz~d^ 


dCt 


f  dcTXz  9(jyZ  9<jz 
\  dx  dy  dz 


+  fz  )  dQ. 


Note  that  because  the  above  equality  holds  for  an  arbitrary  domain  fi,  the  integrands  on 
both  sides  must  be  equal.  Thus,  we  have 


P 


dv ? 


+  vx 


dv? 


+  Vv 


dv? 


+  vz 


dv? 


dt  " z  dx  y  dy  "  ~  dz 


d<xXz  d(JyZ  dazz 

—R - 1 - 7 - 1 - 7T - h  fz, 

dx  dy  dz 


or  written  concisely  as 

Dvz 

P~^  =  V  •  <7.,*  +  fz, 

which  is  the  equation  of  motion  of  a  continuum  in  the  ^-direction.  The  entire  set  for  the 
equations  of  motion  of  a  continuum  in  Eulerian  coordinate  system  is  given  as  follows: 


( dvx  dvx  dvx  dvx 

(  dv.y  dv.y  dVy  <9U y 

p  { w + v‘~d7 + + 

( dv *  dvz  dvz  dvz 


da  XX  dayx  do_  zx  p 

i  “  r  “  i  Jx 


dx 


dy 


dz 


dcrxy  d(Tyy  d  (J  Zy 

- -  +  -77^  +  7  ~  +  fy 


dx 


dy 


dz 


dcrrz  davz  da?z 

■'  +^r1  +  ^1  +  fz 


dx 


dy 


dz 


(2.20) 


We  note  that  (2.20)  is  also  called  Cauchy's  equation  of  motion  or  Cauchy's  momentum 
equation  in  some  literature.  Equation  (2.20)  can  be  written  in  vector  form  as 

P  (jfi  +  (v  •  V)v)  =  V  •  crT  +  f, 

where  cr  is  the  Cauchy  stress  tensor  defined  in  (2.9).  It  is  often  desirable  to  express  these 
equations  of  motion  in  terms  of  displacements  u.  We  find  (because  the  Eulerian  velocity  is 

du 

given  in  terms  of  the  displacement  (2.4)  by  v  =  — ) 


P 


V  •  <tt  +  f . 


(2.21) 


2.2.5  The  Lagrangian  Equations  of  Motion  of  a  Continuum 

Next  we  will  rewrite  (2.20)  in  terms  of  a  Lagrangian  description,  that  is,  we  will  derive  an 
equation  of  motion  in  the  Lagrangian  coordinate  system  (O  —  XYZ  coordinate  system).  Let 
r0  denote  the  boundary  of  Ho  dr  the  initial  (un deformed /material)  configuration,  and  n0  be 
the  outer  normal  vector  on  To-  By  Nanson’s  formula  [40]  we  have 

ndT  =  |H|(Lr1)Tnodr0,  (2.22) 
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where  n0  =  {nox ,  noY ,  noz)T ,  71  is  the  conhguration  gradient  defined  by  (2.5).  Multiplying 
both  sides  of  (2.22)  by  cr1  yields 

<rrndr  =  |A|CTr(A-1)Tn0dr0. 

By  (2.11)  and  the  fact  that  the  Cauchy  stress  tensor  is  symmetric,  we  have 

crTndr  =  Pn0dr0. 


Let  f0  be  the  external  body  force  acting  on  fl0  (fo  =  |^4|f),  let  po(X,  Y,  Z,  t )  be  the  material 
density  in  the  Lagrangian  coordinate  system  (conservation  of  mass  implies  that  p0  —  |7l|p), 
and  V (. X ,  Y,  Z,  t )  be  the  velocity  in  the  Lagrangian  coordinate  system.  Then  we  can  rewrite 
the  resultant  force  in  the  ^-direction  in  the  Eulcrian  coordinate  system  as  the  resultant  force 
in  the  Z  direction  in  the  Lagrangian  coordinate  system,  which  is 


Fqz  —  /  (Pzxnox  +  Pzy^oy  +  Pzznoz)dTQ  +  /  f0ZdVt0. 


■JTo 

Then  by  Gauss’  Theorem  and  the  above  equation  we  find  that 
Fqz  —  ( 


'n0 


'n0 


dPzx  dPzY  dPzz 

'  r\  t  r  ' 


dX 


d  Y 


dZ 


dQ0  +  /  f0Zdn0. 

J  Oq 


We  can  rewrite  Reynolds  transport  theorem  (2.18)  in  the  Lagrangian  coordinate  system  and 
find 


Jt!sda=  Lpitrda  = 


'n  o 


DVz 

Po  dU0. 


DVZ  dVz 

Note  that  — .  =  ^ .  .  Hence,  we  can  rewrite  Newton’s  law  in  the  Lagrangian  coordinate 


system  as 


Dt 


'n  o 


dt 

BVZ 


'n0 


dPZx  ,  dPZY  ,  dPZz  t 

+  7X7  +  +  foZ  dn  0. 


dX 


dY 


dZ 


Note  that  the  above  equality  holds  for  any  Q(>  Thus  we  have 

dVz  dPzx  ,  dPzY  ,  dPzz 


Po- 


dt 


dX 


+ 


dY 


+ 


dZ 


+  foz, 


which  is  the  equation  of  motion  in  the  Z-direction. 
Lagrangian  coordinate  system  are  given  by 

dVx  dPxx  dPXY 
Po  dt  dX  +  dY  + 

dVY  _  dPYX  dPYY 
Po  dt  dX  +  dY  + 

dVz  _  dPzx  dPzY 
Po~df  ~  ~dX~  +  dY  +  ' 


Then  the  equations  of  motion  in  the 


dP 


x  z 


dZ 

9Pyz 

dZ 

dPzz 

dZ 


+  fox 
+  foY 
+  foz, 


(2.23) 
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or,  written  concisely, 


Note  that  V 
given  by 


d\ 

p0~m 


V  •  P  +  f0. 


dU 

~dt' 


Hence,  the  Lagrangian  equations  of  motion  in  terms  of  displacement  is 


<92U 

Po~dt 9 


V  •  P  +  f0. 


(2.24) 


Remark  2.8.  We  note  that  the  equations  of  motion  (2.21)  in  the  Eulerian  (or  moving) 
coordinate  system  are  inherently  nonlinear  independent  of  the  constitutive  law  assumptions 
(discussed  in  the  next  section)  we  might  subsequently  adopt.  On  the  other  hand,  the  La¬ 
grangian  formulation  (2.24)  (relative  to  a  fixed  referential  coordinate  system)  will  yield  a 
linear  system  if  a  linear  constitutive  law  is  assumed.  Thus,  there  are  obvious  advantages 
to  using  the  Lagrangian  formulation  in  linear  theory  (i.e.,  when  a  linear  constitutive  law  is 
assumed) . 


3  Constitutive  Relationships:  Stress  and  Strain 


In  the  preceding  discussions,  we  have  focused  on  relationships  between  displacements  (and 
their  rates)  and  the  stress  tensors.  We  have  also  related  strain  tensors  to  displacements. 
To  complete  our  derivations  of  the  equations  of  motion,  we  must  know  (or  assume)  the 
relationships  (constitutive  “laws”)  between  stress  and  strain.  Constitutive  laws  are  usually 
formulated  based  on  empirical  observations,  and  they  “hold”  for  a  given  material  and  are 
thus  material  dependent.  Moreover,  they  must  be  independent  of  any  referential  coordinate 
system  that  we  choose.  In  addition,  we  must  note  that  the  constitutive  law  describes  an 
ideal  material,  and  it  should  provide  a  close  approximation  to  the  actual  behavior  of  the  real 
material  that  this  constitutive  law  is  intended  to  model. 

The  concept  of  isotropy  is  used  frequently  as  a  simplifying  assumption  in  continuum  me¬ 
chanics,  and  many  useful  materials  are  indeed  isotropic  or  approximately  so.  We  proceed  to 
present  the  formal  definition  of  an  isotropic  tensor  and  isotropic  materials. 

Definition  3.1.  If  a  tensor  has  the  same  array  of  components  when  the  frame  of  reference 
is  rotated  or  reflected  (i.e.,  invariance  under  rotation  or  reflection),  then  it  is  said  to  be  an 
isotropic  tensor.  A  material  whose  constitutive  equation  is  isotropic  is  said  to  be  an  isotropic 
material. 

Remark  3.2.  If  the  tensor  T>ljki  is  isotropic,  then  it  can  be  expressed  in  terms  of  two 
independent  constants  v  and  /i  by 

udjj dki  T  fi (Sjkdji  -f-  &u5jk).  (3.1) 


Since  we  are  interested  in  incorporating  the  constitutive  laws  for  stress  and  strain  into  the 
equations  of  motion,  we  will  only  present  constitutive  laws  for  their  relaxation  forms,  i.e., 
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stress  is  a  function  of  strain  and/or  strain  rate.  The  corresponding  compliance  forms,  i.e., 
the  strain  in  terms  of  stress  and/or  stress  rate,  for  most  of  these  constitutive  laws  can  be 
defined  similarly  by  just  interchanging  the  role  of  stress  and  strain.  For  convenience,  we  will 
suppress  the  spatial  dependence  of  both  stress  and  strain  when  the  constitutive  relationship 
is  given.  Recall  also  that  in  an  infinitesimal  setting  the  stress  tensors  are  all  equivalent;  unless 
noted  otherwise,  we  will  use  a  to  denote  the  stress  in  the  following  discussion  and  assume 
an  infinitesimal  setting.  The  rest  of  this  section  is  outlined  as  follows:  We  first  talk  about 
the  constitutive  equations  used  in  clastic  materials  in  Section  3.1,  and  then  we  present  and 
discuss  a  number  of  constitutive  laws  appearing  in  the  literature  for  the  viscoelastic  materials 
in  Section  3.2. 


3.1  Elastic  Materials 

Elasticity  is  the  physical  property  of  a  material  that  when  it  deforms  under  stress  (e.g., 
external  forces),  but  it  returns  to  its  original  shape  when  the  stress  is  removed.  For  an 
elastic  material,  the  stress-strain  curve  is  the  same  for  the  loading  and  unloading  process, 
and  the  stress  only  depends  on  the  current  strain,  not  on  its  history.  A  familiar  example  of 
an  elastic  material  body  is  a  typical  metal  spring.  Below  we  will  discuss  linear  elasticity  in 
Section  3.1.1  and  then  follow  with  comments  on  nonlinear  elasticity  in  Section  3.1.2. 


3.1.1  Linear  Elasticity 

The  classical  theory  of  elasticity  deals  with  the  mechanical  properties  of  clastic  solids  for 
which  the  stress  is  directly  proportional  to  the  stress  in  small  deformations.  Most  structural 
metals  are  nearly  linear  clastic  under  small  strain  and  follow  a  constitutive  law  based  on 
Hooke’s  law.  Specifically,  a  Hookean  elastic  solid  is  a  solid  that  obeys  Hooke’s  Law 

&ij  Gjkl^kli  (3-2) 

where  Cijki  is  elasticity  tensor.  If  a  material  is  isotropic ,  i.e.,  the  tensor  ctJki  is  isotropic,  then 
by  (3.1)  and  (3.2)  we  have 

(Tij  vSijE kk  T  2 (3.3) 

where  v  and  /i  are  called  Lame’s  parameters.  In  engineering  literature,  the  second  Lame 
parameter  is  further  identified  as  the  shear  modulus. 


3.1.2  Nonlinear  Elasticity 

There  exist  many  cases  in  which  the  material  remains  elastic  everywhere  but  the  stress-strain 
relationship  is  nonlinear.  Examples  are  a  beam  under  simultaneous  lateral  and  end  loads, 
as  well  as  large  deflections  of  a  thin  plate  or  a  thin  shell.  Here  we  will  concentrate  on  the 
hyperelastic  (or  Green  elastic )  material,  which  is  an  ideally  clastic  material  for  which  the 
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strain  energy  density  function  (a  measure  of  the  energy  stored  in  the  material  as  a  result  of 
deformation)  exists.  The  behavior  of  unfilled,  vulcanized  elastomers  often  conforms  closely 
to  the  hyperelastic  ideal. 


Nonlinear  stress-strain  relations  Let  W  denote  the  strain  energy  function,  which  is  a 
scalar  function  of  configuration  gradient  A  defined  by  (2.5).  Then  the  first  Piola-Kirchhoff 
stress  tensor  P  is  given  by 


P 


dW 

~dA' 


or 


Pi 


13 


dW 
dAi:j  ’ 


(3.4) 


where  P^  and  AtJ  are  the  (i,j)  components  of  P  and  A,  respectively.  By  (2.6)  we  can  rewrite 
(3.4)  in  terms  of  Lagrangian  strain  tensor  E. 


P 


A 


dW 

PE  ’ 


or 


Pi 


Aik 


dW 
d  Ekj ' 


(3.5) 


By  (2.13)  and  (3.5)  we  find  that  the  second  Piola-Kirchhoff  stress  tensor  S  is  given  by 


S 


dW 
~dE ’ 


or  Sij 


dW 

dE~j' 


(3.6) 


where  is  the  (i,j)  component  of  S.  By  (2.12)  and  (3.6)  we  find  that  the  Cauchy  stress 
tensor  a  is  give  by 


(7  = 


1  dW 
jA\A~d E 


HT. 


Strain  energy  function  for  isotropic  elastic  materials  For  an  isotropic  material,  the 
configuration  gradient  A  can  be  expressed  uniquely  in  terms  of  the  principal  stretches  (Aj, 
i  =  1,2,3)  or  in  terms  of  the  invariants  (Pi ,  /2,/3)  of  the  left  Cauchy-Green  configuration 
tensor  or  right  Cauchy-Green  configuration  tensor.  Hence,  for  an  isotropic  material,  we  can 
express  the  strain  energy  function  in  terms  of  principal  stretches  or  in  terms  of  invariants. 
Note  that  Ai  =  A2  =  A3  =  1,  I\  =  3,  /2  =  3  and  J3  =  1  in  the  initial  configuration  where  we 
choose  W  =  0.  Hence,  a  general  formula  for  the  strain  energy  function  can  be  expressed  as 

W(Ai,A2,A3) 

OO 

=  a*P  {  [^1(^2  +  ^3)  +  -^2(^3  +  ^1)  +  ^3(^1  +  ^2)]  -  6} 

i,j,k= 0 

or 

OO 

W(h,  I2l  I3)  =  aA1 1  -  3)*(/2  -  3 )j(I3  -  l)k.  (3.8) 

i,j,k= 0 

Then  for  incompressible  materials  (many  rubber  or  elastomeric  materials  are  often  nearly 
incompressible),  |H|  =  1  (which  implies  that  A3A2A3  =  1  and  J3  =  1),  so  (3.7)  can  be  reduced 
to 

OO 

W(Xi,  A2,  A3)  =  ciij  [A^ (A2  +  A3)  +  A^(A^  +  A{)  +  A^(A{  +  A2)  —  6]  (3.9) 

i,j= 0 
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subject  to  A1A2A3  =  1,  and  (3.8)  can  be  reduced  to 


W(h,I2)  =  J2  CiAh  -  Z)V- 2  -  3)J-  (3. tO) 

i,j= 0 

Special  cases  for  (3.10)  include  several  materials: 


•  A  Neo-Hookean  material  for  which  W(/i,/2)  =  cio(/i  —  3), 

•  A  Mooney- Rivilin  (or  Mooney)  material  for  which  W (A,  J2)  =  cw(Ii  —  3)  +  c0i(/2  —  3). 


The  Neo-Hookean  and  Mooney-Rivilin  strain  energy  functions  have  played  an  important 
part  in  the  development  of  nonlinear  elasticity  theory  and  its  application.  The  interested 
reader  should  consult  [23,  40,  44,  53]  and  the  references  therein  for  further  information  on 
hyperclastic  materials. 


3.2  Viscoelastic  Materials 

The  distinction  between  nonlinear  elastic  and  viscoelastic  materials  is  not  always  easily  dis¬ 
cerned  and  definitions  vary.  However  it  is  generally  agreed  that  viscoelasticity  is  the  property 
of  materials  that  exhibit  both  viscous  (dashpot-like)  and  elastic  (spring-like)  characteristics 
when  undergoing  deformation.  Food,  synthetic  polymers,  wood,  soil  and  biological  soft  tis¬ 
sue  as  well  as  metals  at  high  temperature  display  significant  viscoelastic  effects.  Throughout 
this  section,  we  discuss  the  concept  in  a  one-dimensional  formulation,  such  as  that  which 
occurs  in  the  case  of  elongation  of  a  simple  uniform  rod.  In  more  general  deformations  one 
must  use  tensor  analogues  (as  embodied  in  (3.2))  of  the  stress,  the  strain  and  parameters 
such  as  modulus  of  elasticity  and  damping  coefficient. 

In  this  section  we  first  (Section  3.2.1)  introduce  some  important  properties  of  viscoelastic 
materials.  We  then  present  and  discuss  a  number  of  specific  forms  of  constitutive  equa¬ 
tions  proposed  in  the  literature  for  linear  viscoelastic  materials  (Section  3.2.2)  and  those  for 
nonlinear  viscoelastic  materials  (Section  3.2.3). 


3.2.1  Properties  of  Viscoelastic  Materials 

Viscoelastic  materials  are  those  for  which  the  relationship  between  stress  and  strain  depends 
on  time,  and  they  possess  the  following  three  important  properties:  stress  relaxation  (a  step 
constant  strain  results  in  decreasing  stress),  creep  (a  step  constant  stress  results  in  increasing 
strain),  and  hysteresis  (a  stress-strain  phase  lag).  Hysteresis  can  be  seen  from  the  stress- 
strain  curve  which  reveals  that  for  a  viscoelastic  material  the  loading  process  is  different  than 
in  the  unloading  process.  For  example,  the  left  plot  in  Figure  4  illustrates  the  associated 
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stress-strain  curve  for  the  Hookean  elastic  solid,  and  that  in  the  right  plot  of  Figure  4  is 
for  the  Kelvin- Voigt  model  (a  linear  viscoelastic  model  discussed  below  in  Section  3.2.2). 
From  this  figure,  we  see  that  we  can  differentiate  between  the  loading  and  unloading  for 


Figure  4:  Stress  and  strain  curves  during  cyclic  loading-unloading,  (left):  Hookean  elastic 
solid;  (right)  Kelvin- Voigt  material  depicted  by  the  solid  line. 

the  Kelvin- Voigt  material,  but  we  cannot  do  this  for  Hookean  clastic  material.  Thus  the 
Kelvin- Voigt  material  “remembers”  whether  it  is  being  loaded  or  unloaded,  hence  exhibiting 
“hysteresis”  in  the  material. 

Stress  relaxation  In  a  stress  relaxation  test,  a  constant  strain  c0  acts  as  “input”  to  the 
material  from  time  to,  the  resulting  time-dependent  stress  is  decreasing  until  a  plateau  is 
reached  at  some  later  time,  which  is  as  depicted  in  Figure  5.  The  stress  function  G(t) 
resulting  from  the  unit  step  strain  (that  is,  £q  =  1)  is  referred  to  as  relaxation  modulus. 


Figure  5:  Stress  and  strain  histories  in  the  stress  relaxation  test. 


In  a  stress  relaxation  test,  viscoelastic  solids  gradually  relax  and  reach  an  equilibrium  stress 
greater  than  zero  (i.e.,  lim  G(t)  =  G^  >  0),  while  for  viscoelastic  fluids  the  stress  vanishes 

t—>  OO 

to  zero  (i.e.,  lim  G(t)  =  0). 

£—>•00 

Creep  In  a  creep  test,  a  constant  stress  cr0  acts  as  “input”  to  the  material  from  time  t0, 
the  resulting  time-dependent  strain  is  increasing  as  depicted  in  Figure  6.  The  strain  function 
J(t)  resulting  from  the  unit  step  stress  (i.e.,  cr0  =  1)  is  called  creep  compliance. 
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t 


t 


Figure  6:  Stress  and  strain  histories  in  the  creep  test. 


In  a  creep  test,  the  resulting  strain  for  viscoelastic  solids  increases  until  it  reaches  a  nonzero 
equilibrium  value  (i.e.,  lim  J(t)  =  Joo  >  0),  while  for  viscoelastic  fluids  the  resulting  strain 

t—>  OO 

increases  without  bound  as  t  increases. 


Dynamic  loading:  stress-strain  phase  lag  Creep  and  stress  relaxation  tests  are  con¬ 
venient  for  studying  material  response  at  long  times  (minutes  to  days),  but  less  accurate 
at  shorter  times  (seconds  and  less).  Dynamic  tests,  in  which  the  stress  (or  strain)  resulting 
from  a  sinusoidal  strain  (or  stress)  is  measured,  are  often  well-suited  for  filling  out  the  short- 
time  range  of  viscoelastic  response.  We  illustrate  these  ideas  with  a  discussion  of  the  stress 
resulting  from  a  sinusoidal  strain  (as  the  discussion  in  the  strain  resulting  from  sinusoidal 
stress  can  proceed  similarly  by  just  interchanging  the  role  of  stress  and  strain). 

In  a  typical  dynamic  test  carried  out  at  a  constant  temperature,  one  programs  a  loading 
machine  to  prescribe  a  cyclic  history  of  strain  to  a  material  sample  rod  given  by 

e(t)  =  e0sm(ut),  (3.11) 

where  £o  is  the  amplitude,  and  u  is  the  frequency.  The  response  of  stress  as  a  function  of 
time  t  depends  on  the  characteristics  of  the  material  which  can  be  separated  into  several 
categories: 


•  A  purely  clastic  solid. 

For  this  material,  stress  is  proportional  to  the  strain  at  all  the  time,  i.e.,  cr(t)  =  ne(t). 
Hence  with  the  strain  defined  in  (3.11),  the  stress  is  given  by 

a{t)  =  keq  sin(ut). 

We  find  the  stress  amplitude  a0  is  linear  in  the  strain  amplitude  e0;  ao  —  K£o-  The 
response  of  stress  caused  by  strain  is  immediate.  That  is,  the  stress  is  in  phase  with 
the  strain. 


•  A  purely  viscous  material. 

For  this  kind  of  material,  stress  is  proportional  to  the  strain  rate:  a{t) 
strain  defined  in  (3.11),  the  stress  is  then  given  by 


77—.  For  the 
(It 


a(t)  =  rjeooj  cos  (cut). 
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Note  that  cos  (tut) 


sin 


and  thus  we  can  rewrite  the  above  expression  as 


cr(t)  =  r)e0u  sin  {tot  +  ^  . 

The  stress  amplitude  is  linear  in  the  strain  amplitude:  <To  =  y/eocu,  which  is  dependent 
on  the  frequency  cu.  The  stress  is  out  of  phase  with  the  strain,  and  strain  lags  stress 
by  a  90  degree  phase  lag. 


•  A  viscoelastic  solid. 

With  the  sinusoidal  strain  (3.11),  the  stress  as  a  function  of  time  appears  complicated 
in  the  first  few  cycles.  But  a  steady  state  will  eventually  be  reached  in  which  the 
resulting  stress  is  also  sinusoidal,  having  the  same  angular  frequency  cu  but  retarded 
in  phase  by  an  angle  5.  This  is  true  even  if  the  stress  rather  than  the  strain  is  the 
controlled  variable.  The  cyclic  stress  is  written  as 


a(t)  =  <t o  sin  (cut  +  <5),  (3.12) 

where  the  phase  shift  5  is  between  0  and  7t/2,  and  the  stress  amplitude  a0  depends  on 
the  frequency  cu.  By  an  identity  of  trigonometry,  we  can  rewrite  (3.12)  as 

a(t)  =  <r0  cos(h)  sin(c ot)  +  <7 o  sin(h)  cos(cut).  (3.13) 


Thus  the  stress  is  the  sum  of  an  in-phase  response  and  out-of-phase  response. 


Next  we  consider  the  energy  loss  for  a  viscoelastic  material.  Let  l  be  the  length  of  a  rod 
with  cross  sectional  area  a.  When  the  solid  is  strained  sinusoidally,  according  to  (3.11),  the 
solid  elongates  as 

A  l(t)  =  /e0  sin  (cut). 

By  (3.13),  the  force  on  the  rod  is 

F{t)  =  aa{t)  =  acr0  cos (h)  sin(cut)  +  aa0  sin(h)  cos(cut). 

During  a  time  interval  dt,  the  solid  elongates  by  dAl,  and  the  work  done  on  the  rod  is 

F(t)dAl  =  F(t)-jj-dt  =  le0uF{t)  cos(c ot)dt. 

In  one  full  cycle  the  work  done  is 

s'2'jt/lj  /»27t/ci; 

W  =  aleoaoLU  cos(<5)  /  sin(cut)  cos(cut)dt  +  a/£o<AWsin(<5)  /  cos(cut)  cos(cut)dt 

Jo  Jo 

=  Tialeocro  sin(<5). 

Note  that  the  in-phase  components  produce  no  net  work  when  integrated  over  a  cycle,  while 
the  out-of  phase  components  result  in  a  net  dissipation  per  cycle  equal  to: 

W  =  Tialeocro  sin(h). 
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Thus,  for  a  purely  elastic  solid,  the  stress  is  in  phase  with  the  strain  (5  =  0)  and  no  energy 
is  dissipated.  On  the  other  hand,  motion  in  the  viscoelastic  solid  produces  energy  loss. 


It  is  a  common  practice  in  engineering  to  use  complex  variables  to  describe  the  sinusoidal 
response  of  viscoelastic  materials.  Thus,  instead  of  strain  history  (3.11),  we  specify  the 
complex  strain  as  e*  =  £oexP(*^)-  Then  we  obtain  the  following  complex  stress  instead  of 
stress  described  by  (3.12) 

cr*  =  (To  exp(i(ut  +  5)). 

The  above  equation  can  be  rewritten  as 

(7  —  Gr  6  j 


where  G*  is  defined  by 

G*  =  —  exp(i<5)  =  —  cos(<5)  +  i—  sin(<5).  (3-14) 

So  So  So 

The  characteristic  parameter  G*  is  referred  to  as  the  complex  dynamic  modulus.  We  denote 
the  real  part  of  G*  by  G'  and  the  imaginary  part  of  G*  by  G" .  That  is, 


G*  =  G'  +  iG",  (3.15) 

where  G'  =  —  cos(<5)  and  G"  =  —  sin(<5).  The  coefficient  G'  is  called  the  storage  modulus 
£o  eo 

(a  measure  of  energy  stored  and  recovered  per  cycle)  which  corresponds  to  the  in-phase 
response,  and  G"  is  the  loss  modulus  (a  characterization  of  the  energy  dissipated  in  the 
material  by  internal  damping)  corresponding  to  the  out-of  phase  response.  The  in-phase 
stress  and  strain  results  in  elastic  energy,  which  is  completely  recoverable.  The  7t/2  ont-of- 
phase  stress  and  strain  results  in  the  dissipated  energy. 


Remark  3.3.  The  relationship  between  the  two  transient  functions,  relaxation  modulus 
G(t)  and  creep  compliance  Jit),  for  a  viscoelastic  material  is  given  by 

J(s)G(t  —  s)ds  =  t. 

The  relationship  between  the  relaxation  modulus  G{t)  and  dynamic  modulus  functions  G ' 
and  G"  is  given  by 


G» 


=  Goo  +  u  (G(t)  —  Goo)  sin (ut)dt, 


=  uj  ( G(t )  —  Goo)  cos (ut)dt. 


Jo 

The  relations  may  be  inverted  to  obtain 


G(t)  —  Goo  H — 


G'(u)  -  Gr 


7T 


C 0 


■  sin(ut)dco, 


or 


G[t)  —  Goo  H — 


G"{u) 


7 r 


u J 


cos(ujt)du. 


The  interested  reader  can  refer  to  the  recent  text  [31]  for  further  information. 
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3.2.2  Linear  Viscoelasticity 


The  characteristic  feature  of  linear  viscoelastic  materials  is  that  the  stress  is  linearly  pro¬ 
portional  to  the  strain  history,  and  it  is  important  to  note  that  the  property  of  linearity  of 
response  does  not  refer  to  the  shape  of  any  material  response  curve.  Linear  viscoelasticity  is 
usually  applicable  only  for  small  deformations  and/or  linear  materials.  Thus,  infinitesimal 
strain  theory  should  be  employed  for  this  case.  There  are  two  standard  approaches  that  have 
been  used  to  develop  constitutive  equations  for  the  linear  viscoelastic  materials:  mechanical 
analogs  and  the  Boltzmann  superposition  principle. 


Mechanical  analogs  Linear  viscoelastic  behavior  can  be  conceived  as  a  linear  combina¬ 
tions  of  springs  (the  clastic  component)  and  dashpots  (the  viscous  component)  as  depicted 
by  Figure  7.  The  elastic  component  is  described  by 


K  7] 


Figure  7:  (left):  Schematic  representation  of  the  Hookean  spring;  (right)  Schematic  repre¬ 
sentation  of  Newtonian  dashpot. 


<7  =  K£, 


de  1  da 
dt  k  dt 

where  a  is  the  stress,  £  is  the  strain  that  occurs  under  the  given  stress,  and  k  is  the  elastic 
modulus  of  the  material  with  units  N/m2 .  The  viscous  component  is  modeled  by 


de 

a  =  'dt' 


where  r/  is  the  viscosity  of  the  material  with  units  N  ■  s/m2.  The  mechanical  analogues 
approach  results  in  a  linear  ordinary  differential  equation  with  constant  coefficients  relating 
stress  and  its  rates  of  finite  order  with  strain  and  its  rates  of  the  form 


da  d2a 

a  +  aidt+a27¥ 


dna  de  d2e  dne 

an—r—  =  b0e  +  +b\  —  +  o-2— r-r  +  . . .  +  bn— 


dtr 


dt 


'  dt 2 


'  dtri 


(3.16) 


In  (3.16)  the  constant  coefficients  a*,  bt,  i  —  0, 1,  2, ...  ,n  are  related  to  the  elastic  modulus 
and  viscosity  of  the  material  which  are  usually  determined  from  physical  experiments.  A 
complete  statement  of  the  constitutive  equation  obtained  from  using  of  mechanical  analogs 
then  consists  of  both  an  equation  of  the  form  (3.16)  and  a  set  of  appropriate  initial  conditions. 
In  addition,  we  see  that  it  is  not  convenient  to  directly  incorporate  this  general  model  (3.16) 
into  the  equation  of  motion. 
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The  three  basic  models  that  are  typically  used  to  model  linear  viscoelastic  materials  are  the 
Maxwell  model,  the  Kelvin- Voigt  model  and  the  standard  linear  solid  model.  Each  of  these 
models  differs  in  the  arrangement  of  these  “springs”  and  “dashpots” . 


The  Maxwell  model  The  Maxwell  model  is  represented  by  a  purely  viscous  damper  and 
a  purely  clastic  spring  connected  in  series  as  depicted  in  Figure  8.  Because  the  spring  and 

at  ^ 

- mm — 3 — •<* 


Figure  8:  Schematic  representation  of  the  Maxwell  model. 


the  dashpot  are  subject  to  the  same  stress,  the  model  is  also  known  as  an  iso-stress  model. 
The  total  strain  rate  is  the  sum  of  the  clastic  and  the  viscous  strain  contributions,  so  that 


a  1  da 
rj  k  dt 


de 

dt 


(3.17) 


We  consider  the  stress  relaxation  function  and  creep  function  for  the  Maxwell  model  (3.17). 
The  stress  relaxation  function  corresponds  to  the  relaxation  that  occurs  under  an  imposed 
constant  strain,  given  eft)  =  £0 Hft  —  to)  and  cr(0)  =  0,  where  Hft)  is  the  Heaviside  step 
function  (also  called  the  unit  step  function  in  the  literature),  t0  >  0.  The  solution  aft)  to 
(3.17)  is  the  relaxation  function.  With  this  strain  function,  (3.17)  can  be  written  as 

a  1  da  r/  . 

— I - 7 7  —  £o  oft  —  t0), 

77  k  dt 

where  5  is  the  Dirac  delta  function.  Let  a(s)  =  Jf{a(t)}(s)  where  Jzf  denotes  the  Laplace 
transform.  Then  taking  the  Laplace  transform  of  both  sides  of  the  above  differential  equation 
we  obtain 

g  *0$  g  tos 

°(s)  =  £ot~t;  = 

r\  '  k  r) 

Taking  the  inverse  Laplace  transform  one  finds  that 

aft)  =  exp  -  t0)^j  £0 Hft  -  to). 

The  stress  relaxation  function  for  the  Maxwell  model  (3.17)  is  illustrated  in  Figure  9  (compare 
with  Figure  5). 

The  creep  function  corresponds  to  the  creep  that  occurs  under  the  imposition  of  a  constant 
stress  given  by  aft)  =  a0H(t  —  t0)  and  e(0)  =  0;  the  solution  eft)  to  (3.17)  is  the  creep 
function.  With  this  stress  function,  (3.17)  can  be  written  as 

^  =  —Hft  -  to)  +  —  6 ft  -  t0). 
dt  rj  ft 
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The  Kelvin-Voigt  model  The  Kelvin- Voigt  model,  also  known  as  the  Voigt  model,  con¬ 
sists  of  a  Newtonian  damper  and  a  Hookean  elastic  spring  connected  in  parallel,  as  depicted 
in  Figure  11.  Because  the  two  elements  are  subject  to  the  same  strain,  the  model  is  also 


Figure  11:  Schematic  representation  of  the  Kelvin-Voigt  model. 


known  as  an  iso-strain  model.  The  total  stress  is  the  sum  of  the  stress 
stress  in  the  dashpot,  so  that 


de 

a  =  H£  +  77  —  . 

dt 


in  the  spring  and  the 
(3.18) 


We  first  consider  the  stress  relaxation  function  and  creep  function  for  the  Kelvin-Voigt 
model  (3.18).  The  stress  relaxation  function  corresponds  to  the  solution  of  (3.18)  when 
e(t )  =  e0H(t  —  t0)  and  cr(0)  =  0.  We  find 

cr(t)  =  K£0H(t  -  t0)  +  rj£0S(t  -  t0)- 

This  stress  relaxation  function  for  the  Kelvin-Voigt  model  (3.18)  is  illustrated  in  Figure  12. 


<T 

- 

dt) 

0 

f0  1  L  * 

Figure  12:  Stress  relaxation  function  for  the  Kelvin-Voigt  model. 


The  creep  function  again  is  the  solution  £{t)  to  (3.18)  corresponding  to  cr(t)  =  o-0H(t  —  to) 
and  e(0)  =  0.  We  find 

d£  . 

K£  +  T,~di  =  ~  to 

Let  £(s)  =  2z?{e(t)}(,s),  or  in  terms  of  the  Laplace  transform  we  have 

g  tos  g  tos 

s  s  +  n/rj 


D  —  CQ  s 


£(s)  =  a  o 


_ _  °o 

s(k  +  rjs)  k 
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Using  the  inverse  Laplace  transform  we  obtain 


e{t) 


1 

K 


1  —  exp 


a0H(t  - 10). 


The  corresponding  creep  function  for  the  Kelvin- Voigt  model  (3.18)  is  illustrated  in  Figure 
13  (compare  with  Figure  6). 


Figure  13:  Creep  function  of  Kelvin- Voigt  model. 


Thus  we  find  that  the  Kelvin- Voigt  model  is  extremely  accurate  in  modelling  creep  in  many 
materials.  However,  the  model  has  limitations  in  its  ability  to  describe  the  commonly  ob¬ 
served  relaxation  of  stress  in  numerous  strained  viscoelastic  materials. 


Finally  we  may  obtain  the  storage  modulus  G'  and  loss  modulus  G"  for  the  Kelvin- Voigt 
model  (3.18).  Letting  e(t)  =  £0exp(itot)  and  a(t)  =  a0  exp(i(tot  +  5))  and  substituting  £  and 
a  into  (3.18),  we  find 

G*  =  —  exp(i<5)  =  k  +  irjuj. 
eo 


Hence,  by  (3.15)  we  have 


G'  =  K,  G"  =  T]U. 


Remark  3.4.  Because  of  one  of  our  motivating  applications  (discussed  in  Section  4  below), 

we  are  particularly  interested  in  the  elastic/viscoelastic  properties  of  soil.  Hardin  in  [29] 

presented  an  analytical  study  of  the  application  of  the  Kelvin- Voigt  model  to  represent  dry 

soils  for  comparison  with  test  results.  From  this  study  he  found  that  the  Kelvin- Voigt  model 

satisfactorily  represented  the  behavior  of  sands  in  these  small- amplitude  vibration  tests  if  the 

viscosity  r)  in  the  model  was  treated  as  varying  inversely  with  the  frequency  to  to  maintain 
•quo 

the  ratio  —  constant.  Since  Hardin’s  work,  describing  soils  as  a  Kelvin- Voigt  material  has 

K 

become  accepted  as  one  of  the  best  ways  in  soil  dynamics  of  calculating  wave  propagation 
and  energy  dissipation.  The  Kelvin- Voigt  model  also  governs  the  analysis  of  standard  soil 
tests,  including  consolidation  and  resonant-column  tests.  The  interested  reader  can  refer  to 
[36,  37,  43]  as  well  as  the  references  therein  for  more  information  of  the  model’s  historical 
backgrounds  and  examples  of  its  application  in  soil  dynamics. 


It  is  interesting  to  note  that  the  author  in  [13]  showed  that  the  Kelvin- Voigt  model  can  be 
used  to  describe  the  dynamic  response  of  the  saturated  poroelastic  materials  that  obey  the 
Biot  theory  (the  two-phase  formulation  of  Biot  in  [15]).  This  viscoelastic  model  is  simpler  to 
use  than  poroelastic  models  but  yields  similar  results  for  a  wide  range  of  soils  and  dynamic 
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loadings.  In  addition,  the  author  in  [36]  developed  a  model,  the  Kelvin- Voigt-Maxwell-Biot 
model,  that  splits  the  soil  into  two  components  (pore  fluid  and  solid  frame)  where  these  two 
masses  are  connected  by  a  dashpot  which  can  then  be  related  to  permeability.  In  addition, 
a  mapping  between  the  Kelvin- Voigt  model  and  the  Kelvin- Voigt-Maxwell-Biot  model  is 
developed  in  [36]  so  that  one  may  continue  to  use  the  Kelvin- Voigt  model  for  saturated  soil. 


The  standard  linear  solid  (SLS)  model  The  standard  linear  solid  model,  also  known 
as  the  Kelvin  model  or  three- element  model ,  combines  the  Maxwell  Model  and  a  Hookean 
spring  in  parallel  as  depicted  in  Figure  14. 


— wwr — 1 


Figure  14:  Schematic  representation  of  the  standard  linear  model. 


The  stress-strain  relationship  is  given  by 

da  de 

a  +  T‘Tt=  Kr(£  +  >' 

/ j |  Kr  T  K\  ,  i  •  i 

where  t£  =  — ,  and  Ta  =  r)i -  (from  which  can  be  observed  that  rCT  >  t£ 

KrK  1 


(3.19) 


The  stress  relaxation  function  and  the  creep  function  for  the  standard  linear  model  (3.19) 
are  obtained  in  the  usual  manner.  As  usual,  the  stress  relaxation  function  a{t)  is  obtained 
by  solving  (3.19)  with  e{t)  =  £0 H(t  —  to)  and  cr(0)  =  0.  We  find 


Cr(t)  +  T£^-  =  Kr£0H(t  -  t0)  +  KrTa£0S(t  -  t0), 
dt 


or  in  terms  of  the  Laplace  transform 


b(s) 


0  —  tQS 


Kr£o 


s(l  +  t£s) 


g  tos 

K,r7-(j£ o~  j 

1  +  r£s 


g  ^0  5 

0 -  0 

S 


g  tos 
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Thus  we  find 


cr(t)  =  Kr£0H(t  -  t0)  +  nr£0  (  —  -  1  exp  (-(t  -  t0)/rs)  H(t  -  t0) 


rf 


=  Kr 


Trr 


Tr 


1  +  —  -  1  exp  (-(t  -  t0)/r£) 


£oH(t  —  to) 


=  [nr  +  ki  exp  {-(t  -  to)/ t£)\  e0H(t  -  t0). 

This  stress  relaxation  function  of  standard  linear  model  (3.19)  is  illustrated  in  Figure  15. 


Figure  15:  Stress  relaxation  function  of  for  the  standard  linear  model. 


The  creep  function  is  the  solution  of  (3.19)  for  e(t)  given  a(t)  =  cr0H(t  —  t0)  and  e(0) 
Using  the  same  arguments  as  above  in  finding  the  stress  function,  we  have 


£{t)  = 


Hr 


1  + 


Tr 


Trr 


1  exp  (-(t  -  t0) /T^) 


CT0H(t  -  t0). 


The  creep  function  of  standard  linear  model  (3.19)  is  illustrated  in  Figure  16. 


0. 


Figure  16:  Creep  function  for  the  standard  linear  model. 


We  therefore  see  that  the  standard  linear  model  is  accurate  in  predicating  both  creep  and 
relaxation  responses  for  many  materials  of  interest. 


Finally  the  usual  arguments  for  the  standard  linear  model  lead  to 

(To  r.rr  «V(1  +  TaT£UJ2)  .KriTry  -T£)UJ 
G  =  — exp  (id)  =  - - v^- 

£0  1  +  (Teu)2  1  +  (TeUj)2 

By  definition  of  the  storage  and  loss  modulus,  we  thus  have 


G'  = 


Kr{l  +  TaTEUJ2) 

1  +  (t€u)2 


G"  = 


Kr(Ta  ~  Te)0J 

1  +  (t€u)2 
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Remark  3.5.  It  was  demonstrated  in  [41]  that  the  relaxation  behavior  of  compressed  wood 
can  be  adequately  described  by  the  standard-linear-model. 


The  Boltzmann  superposition  model  A  general  approach  widely  used  to  model  linear 
viscoelastic  materials  is  due  to  Boltzmann  (1844-1906)  and  is  called  the  Boltzmann  superpo¬ 
sition  model  or  simply  the  Boltzmann  model.  If  the  origin  for  time  is  taken  at  the  beginning 
of  motion  and  loading  (i.e. ,  cr(0)  =  0  and  e(0)  =  0),  then  the  stress-strain  law  is  given  by 

a(t)  —  Kre(t)  +  f  Kit  —  s)—  7  - ds,  (3.20) 

Jo  ds 

where  Kr  represents  an  instantaneous  relaxation  modulus,  and  K  is  the  “gradual”  relaxation 
modulus  function.  The  relaxation  modulus  function  G(t)  for  the  Boltzmann  model  (3.20) 
is  given  by  G(t)  =  ftr  +  K(t).  We  find  that  when  model  (3.20)  is  incorporated  into  force 
balance  laws  (equation  of  motion),  it  results  in  integro-partial  differential  equations  which 
are  most  often  phenomenological  in  nature  as  well  as  being  computationally  challenging  both 
in  simulation  and  control  design. 

Remark  3.6.  Any  constitutive  equation  of  form  (3.16),  along  with  the  appropriate  initial 
conditions,  can  be  expressed  in  either  the  form  (3.20)  or  its  corresponding  compliance  form 
(i.e.,  the  strain  in  terms  of  the  stress  and  stress  rate).  On  the  other  hand,  a  constitutive 
equation  of  form  (3.20)  can  be  reduced  to  the  form  (3.16)  if  and  only  if  the  stress  relaxation 
modulus  satisfy  specific  conditions.  A  detailed  discussion  of  this  can  be  found  in  [27].  For 
example,  the  Maxwell,  Kelvin- Voigt,  and  standard  linear  models  can  be  expressed  in  the 
Boltzmann  formulation.  The  choice  of  parameters  in  (3.20)  to  yield  these  models  are  readily 
verified  to  be: 


V 

•  The  Maxwell  model:  nr  =  0,  and  Kit  —  s)  =  ft  exp  (—it  —  s)  It)  with  r  =  — ; 

ft 

•  The  Kelvin- Voigt  model:  ftr  =  ft,  and  K(t  —  s)  =  rj5(t  —  s); 


Trr 


The  standard  linear  model:  K(t  —  s)  —  —  1 J  ftrexp  (— (t  —  s)/re ),  which  can  be 

Tji 

simplified  as  Kit  —  s)  =  fti  exp  (  —  (t  —  s)/re)  with  t£  =  — . 

fti 


Internal  variable  approach  Other  special  cases  of  (3.20)  that  are  often  encountered 
include  a  generalization  of  the  single  spring-dashpot  paradigm  of  Figure  14  (the  standard 
linear  model)  to  one  with  multiple  spring-dashpot  systems  in  parallel  as  depicted  in  Figure 
17.  This  is  a  generalized  standard  linear  model  (also  referred  to  as  generalized  Maxwell 
model  or  Wiechert  model  or  Maxwell-  Wiechert  model  in  the  literature)  that  results  in  the 
relaxation  response  kernel 

n 

Kit  -s)  =  Y^  Kj  exp  (-(t  -  s)/Tj) ,  (3.21) 

3= 1 
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Figure  17: 


Schematic  representation  of  the  generalized  standard  linear  model. 


where  t?-  =  — ,  j  =  1,2 , ,n.  We  observe  that  the  Boltzmann  formulation  with  relaxation 

Kj 

response  kernel  (3.21)  is  equivalent  to  the  formulation 


a(t)  =  Kre{t)  +  E  KAW,  (3-22) 

3= 1 


where  the  satisfy  the  following  ordinary  differential  equations 


d£j(f) 

dt 


de{t ) 
dt  ’ 


e(°)  =  0,  j 


1,2,.. 


n. 


(3.23) 


The  formulation  (3.22)  with  (3.23)  is  sometimes  referred  to  as  an  internal  variable  model  (e.g. 
see  [7,  12])  because  the  variables  e3  can  be  thought  of  as  “internal  strains”  that  are  driven  by 
the  instantaneous  strain  according  to  (3.23)  and  that  contribute  to  the  total  stress  via  (3.22). 
Hence,  we  can  see  that  this  internal  variable  modeling  leads  to  an  efficient  computational 
alternative  for  the  corresponding  integro-partial  differential  equation  models  involving  (3.20). 
In  addition,  we  note  that  this  approach  provides  a  “molecular”  basis  for  the  models  as  it  can 
be  thought  of  as  a  model  for  a  heterogeneous  material  containing  multiple  types  of  molecules 
[7,  12,  18,  30],  each  possessing  a  distinct  relaxation  property  characterized  by  a  relaxation 
parameter  t3.  For  a  comparison  of  models  of  viscoelastic  damping  via  hysteretic  integrals 
versus  internal  variable  representations,  see  [7]  and  the  references  therein. 


Moreover,  we  see  that  model  (3.22)  with  (3.23)  can  be  readily  viewed  as  a  special  case  of 
or  even  an  approximation  to  models  with  a  continuum  of  relaxation  times  (see  [12]  and 
the  references  therein).  These  models  have  proved  useful  in  a  wide  range  of  viscoelastic 
materials.  The  corresponding  stress-strain  laws  have  the  form 

a(t]  V)  =  Kre{t)  +  +  7  J  £i (t;  r)d'P(r),  (3.24) 

where  V  is  a  probability  distribution  over  the  set  T  of  possible  relaxation  parameters  r,  and 
£i  (i;  r)  satisfies,  for  each  r  G  T, 


efei(i;  r) 
dt 


+  -ei(t;  r) 
r 


de{t) 

dt 


£i(0;t)  =  0. 


(3.25) 
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The  approach  embodied  in  equations  (3.24)-(3.25)  offers  a  computationally  tractable  alter¬ 
native  for  linear  materials  with  a  continuum  of  relaxation  times. 

Remark  3.7.  The  generalized  standard  linear  model  (i.e.,  model  (3.22)  combined  with 
(3.23))  with  different  n  has  been  successfully  used  to  describe  the  stress  relaxation  behavior 
of  a  variety  of  foods.  For  example,  in  [51]  this  model  with  n  —  2  was  successfully  used  to 
describe  the  stress  relaxation  of  lipids  such  as  beeswax,  candelilla  wax,  carnauba  wax  and 
a  high-melting  milkfat  fraction.  A  comprehensive  study  on  the  ability  of  the  generalized 
Maxwell  model  to  describe  the  stress  relaxation  behavior  of  solid  food  is  presented  in  [17]. 
In  this  study,  five  different  food  matrices  (agar  gel,  meat,  ripened  cheese,  “mozzarella”  cheese 
and  white  pan  bread)  were  chosen  as  representatives  of  a  wide  range  of  foods,  and  results 
verify  that  the  proposed  model  satisfactorily  fits  the  experimental  data. 


3.2.3  Nonlinear  Viscoelasticity 

For  many  materials  the  linear  models  given  in  Section  3.2.2  are  inadequate  to  describe 
experimental  data.  In  particular,  shape  memory  alloys  such  as  Nitinol  (a  nickel-titanium) 
and  CuZnAl  (a  cooper  zinc  aluminum  aloy),  biological  soft  tissue  and  highly  filled  rubber 
exhibit  significant  nonlinear  hysteric  behavior  such  as  that  depicted  by  the  experimental  data 
for  highly  filled  rubber  given  in  Figure  18.  Nonlinear  viscoelastic  behavior  is  usually  exhibited 
when  the  deformation  is  large  or  if  the  material  changes  its  properties  under  deformations. 


Figure  18:  Experimental  stress-strain  curves  for  (1)  unfilled,  (2)  lightly  filled  and  (3)  highly 
filled  rubber  in  tensile  deformations. 


The  theory  of  nonlinear  viscoelasticity  has  attracted  the  attention  of  a  large  number  of 
investigators  over  the  past  century  (e.g.  [18,  19,  22,  24,  33,  38,  47,  52]).  Two  types  of 
models  for  stress-strain  relationships  can  be  found  in  the  literature.  One  is  based  on  the 
phenomenological  mechanical  behavior  of  the  materials  (that  is,  the  form  of  constitutive 
equations  is  not  based  on  the  explanation  of  how  these  properties  arise  from  the  underlying 
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microscopic  structure).  For  example,  Green  and  Rivlin  in  [26]  constructed  a  multiple  integral 
constitutive  equation,  which  is  arranged  as  a  series  in  which  the  nth  term  is  of  degree  n  in 
the  strain  components.  A  multiple  integral  constitutive  equation,  arranged  in  a  series,  was 
also  developed  by  Pipkin  and  Rogers  in  [42],  in  which  the  first  term  gives  the  results  of  a 
one-step  test  (the  stress  due  to  a  step  change  of  strain),  and  whose  nth  term  represents 
a  correction  due  to  the  nth  step.  The  interested  reader  can  refer  to  [19,  55]  for  recent 
historical  overviews  on  these  phenomenological  models  as  well  as  the  mathematical  issues 
underlying  the  formulations  of  these  models.  The  other  type  of  model  entails  formulations 
based  on  the  molecular  mechanisms  underlying  the  response.  For  example,  in  [18]  Doi 
and  Edwards  developed  a  “reptation”  model  for  concentrated  solutions  and  polymer  melts 
which  is  based  on  the  assumption  that  an  entangled  polymer  molecule,  the  chain,  slides 
(or  reptates)  through  a  “tube”  whose  contours  are  defined  by  the  locus  of  entanglements 
with  neighboring  molecules.  The  polymer  chain  is  free  to  diffuse  along  the  tube  axis  but 
cannot  move  perpendicularly  to  the  tube  as  other  molecules  restrict  such  movement.  In  [14] 
the  polymer  melt  is  modeled  as  a  collection  of  interacting  bead-rod  chains.  This  model, 
referred  to  as  the  Curtiss-Bird  model,  is  similar  to  the  Doi-Edwards  model,  but  it  is  based 
on  a  systematic  kinetic  development  and  does  not  use  the  phenomenological  constraints  of 
a  chain  in  a  tube.  In  addition,  the  comparison  results  of  the  Curtiss-Bird  model  and  the 
Doi-Edwards  model  in  [14]  show  that  the  the  Curtiss-Bird  model  is  more  accurate  than 
Doi-Edwards  model  in  predicting  the  nonlinear  behavior.  A  review  of  these  molecular  types 
of  models  as  well  as  the  relationship  between  the  nonlinear  viscoelasticity  and  molecular 
structure  is  given  in  [39]. 

In  the  following  sections  we  restrict  our  discussions  to  those  models  that  have  been  employed 
and  developed  by  our  group  for  understanding  the  dynamic  response  of  the  highly  filled 
rubber  and  propagation  of  arterial  stenosis  induced  shear  waves  in  composite  biotissue.  The 
forms  of  these  constitutive  equations  can  be  easily  incorporated  into  the  equation  of  motion 
so  that  we  can  numerically  solve  a  partial  differential  equation  with  appropriate  boundary 
and  initial  conditions  to  understand  the  dynamic  behavior  of  materials. 


Modified  Boltzmann  superposition  model  The  most  direct  formulation  to  treat  non¬ 
linear  viscoelasticity  is  one  based  on  generalizing  the  Boltzmann  superposition  model  (3.20) 
to  a  corresponding  nonlinear  version.  That  is,  one  allows  nonlinear  instantaneous  strain  as 
well  as  nonlinear  strain  rate  dependence.  This  “modified”  superposition  principle  was  first 
suggested  in  [32]  where  it  was  observed  the  creep  behavior  of  fibres  could  be  separated  into 
time  and  stress-dependent  parts,  and  then  considered  in  [21]  for  a  multiple-step  test.  One 
form  of  this  approach  has  been  employed  in  some  of  our  earlier  efforts  [9,  10]  for  modeling 
hysteretic  damping  in  elastomers  and  is  given  by 

r1  d, 

cr(t)  =  ge{e(t))  +  cDe(t)  +  J  K(t  -  s)—gv(e(s),e(s))ds,  (3.26) 

where  £  is  the  infinitesimal  strain,  K  is  the  convolution  memory  kernel,  and  ge  and  gv  are 
nonlinear  functions  accounting  for  the  clastic  and  viscoelastic  responses  of  the  materials, 
respectively.  As  explained  in  [9],  our  nonlinear  materials  undergoing  large  deformations 
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required  the  use  of  finite  (as  opposed  to  infinitesimal)  strain  theories.  However,  since  the 
nonlinearity  between  the  stress  and  finite  strain  is  an  unknown  and  is  to  be  estimated 
(using  inverse  problem  algorithms)  and  since  the  finite  strain  can  be  expressed  in  terms  of 
known  nonlinearities  as  a  function  of  the  infinitesimal  strain  (at  least  in  the  problems  of 
interest  here),  one  can  effectively  formulate  the  problem  as  one  of  estimating  the  unknown 
nonlinearity  between  stress  and  infinitesimal  strain.  Thus  one  can  develop  nonlinear  models 
for  large  deformation-related  stress  in  terms  of  infinitesimal  strain. 

Again,  this  form  of  model  (3.26),  when  incorporated  into  force  balance  laws  (equation  of 
motion),  results  in  integro-partial  differential  equations  which  are  computationally  challeng¬ 
ing  both  in  simulation  and  control  design.  This  motives  us  to  employ  the  internal  variable 
approach  described  below  as  an  alternative  for  computation. 


Fung’s  quasi-linear  viscoelastic  formulation  This  formulation  can  be  found  in  Fung 

[24]  in  the  context  of  efforts  to  model  soft  tissue,  which  is  often  characterized  by  a  constant 
hysteresis  over  a  wider  frequency  range.  In  this  formulation,  a  continuous  spectrum  over  a 
finite  range  (ti,T2)  is  used  to  describe  this  phenomenon  instead  of  using  a  finite  spectrum 
{ti,T2,  . . . ,  t„}  associated  with  relaxation  kernels  of  the  form  given  in  (3.21).  Fung’s  model 
is  given  by 


.  .  .  dcre(X(r)) 

a(t)  —  /  K(t  —  t) - - - -dr. 

./n  dT 


(3.27) 


Here  the  stress  a  is  the  force  in  the  deformed  configuration  divided  by  the  cross-sectional  area 
of  the  speciman  at  the  zero  stress  state  (i.e.,  the  scalar  version  of  the  first  Piola-Kirchhoff 
stress  tensor),  A  is  the  ratio  of  the  length  of  the  specimen  stretched  under  the  load  divided 
by  the  initial  length  at  the  zero  stress  state,  ae  is  the  “elastic”  response  to  this  elongation, 
and  K(t)  is  the  reduced  relaxation  function  (“reduced”  refers  to  that  K(t)  is  normalized 
such  that  A' (0)  =  1).  Fung  proposes  a  kernel  of  the  form 


K{t) 


1  +  cJ^expj-t/r)  dr 
1  +  c1ii(t2/ti) 


(3.28) 


where  c  represents  the  degree  to  which  viscous  effects  are  present,  and  T\  and  r2  represent 
fast  and  slow  viscous  time  phenomena. 


Since  its  introduction,  this  quasi-linear  viscoelastic  (QLV)  theory  of  Fung  has  been  applied 
successfully  in  stress-strain  experiments  to  several  types  of  biological  tissue.  A  benefit  to 
using  (3.27)  as  a  constitutive  equation  is  that,  unlike  simpler  models  for  viscoelasticity, 
it  allows  for  the  consideration  of  a  continuous  spectrum  (e.g.,  see  the  discussions  in  [24]) 
of  relaxation  times;  this  is  also  true  of  the  probabilistic-based  internal  variable  approach 
developed  in  [8]  and  described  below.  The  need  for  a  continuum  of  relaxation  times  in  certain 
materials  was  observed  many  years  ago  (e.g.,  see  [20,  28,  49]).  While  Fung’s  theory  has  been 
successfully  employed  for  fitting  hysteretic  stress-strain  curves,  for  control  applications  one 
is  interested  in  using  it  in  a  full  dynamical  model.  Unfortunately,  the  QLV,  as  presented  by 
Fung,  leads  to  exceedingly  difficult  computations  within  full  dynamical  partial  differential 
equations,  especially  in  estimation  and  control  problems. 
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Internal  variable  approach  To  overcome  the  computational  challenges  present  in  the 
modified  Boltzmann  superposition  model  and  Fung’s  quasi- linear  viscoelastic  model,  we  turn 
to  an  internal  variable  approach  in  an  alternative  way  to  formulate  the  constitutive  equation. 


We  observe  that  if  the  relaxation  response  kernel  function  K  in  model  (3.26)  is  defined  as  in 
(3.21),  then  (3.26)  can  be  written  in  terms  of  internal  strains  similar  to  (3.22)  with  (3.23): 


cr(f)  =  ge{e{t))  +  cDe(t)  +  ^  Kj£j{t),  (3.29) 

3= 1 


where  the  internal  strains  £3  satisfies 

+  —  £j(t)  =  ffv(e(t),£(t)),  £j(0)  =  0,  j  =  1,2, .  .  .,n.  (3.30) 

at  Tj 

Through  comparison  with  experimental  data,  in  [10]  it  was  discovered  that  the  best  fit  to 
filled  elastomer  data  occurs  when  ge  and  gv  are  cubic.  Moreover,  model  (3.29)  with  (3.30) 
can  be  readily  generalized  to  models  with  a  continuum  of  relaxation  times  (see  [12]  and  the 
references  therein).  These  models  have  been  useful  in  a  wide  range  of  viscoelastic  materials. 
The  corresponding  stress-strain  laws  have  the  form 


cr(t;V)  =  ge(e(t))  +  cDi(t)  £i (t',r)dV(r),  (3.31) 

where  V  is  a  probability  distribution  over  the  set  T  of  possible  relaxation  parameters  r,  and 
£\ (£;  r)  satishes,  for  each  r  G  T, 

+  -£l(t-T)  =  gv(£(t),£(t)),  £i(0;t)  =  0.  (3.32) 

Cit  T 


We  also  observe  that  if  the  reduced  relaxation  function  K(t)  in  Fung’s  QLV  model  (3.27)  is 

n 

defined  as  in  (3.21)  with  k3  =  1,  then  model  (3.27)  is  equivalent  to  the  following  internal 

k= 1 

strain  variable  formulation 

n 

a(t)  =  '52Kj£j(t),  (3.33) 

3=  1 

with  the  internal  strains  Ej  satisfying 


d£j  (■ t ) 
dt 


dae(\(t)) 

dt 


£j(0)  =  0,  j 


1,2, ...  ,n. 


(3.34) 


Various  internal  strain  variable  models  are  investigated  in  [1]  and  a  good  agreement  is  demon¬ 
strated  between  a  two  internal  strain  variable  model  (i.e.  n  =  2)  and  undamped  simulated 
data  based  on  the  Fung  kernel  (3.28).  The  corresponding  stress-strain  laws  for  a  generaliza¬ 
tion  of  these  models  to  have  a  continuum  of  relaxation  times  have  the  form 


*(t;  V)  =  J^£i(t;r)dV(r), 


(3.35) 
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where  once  again  V  is  a  probability  distribution  over  the  set  T  of  possible  relaxation  param¬ 
eters  r,  and  £\  (t;  r)  satishes,  for  each  r  G  T, 


d£i(t]  r ) 
dt 


+  -ei(t;r) 
r 


dae(X(t)) 

dt 


£i(0;t)  =  0. 


(3.36) 


Molecular  type  of  constitutive  equations  Even  though  model  (3.29)  combined  with 
(3.30)  provides  a  reasonable  fit  to  the  experimental  data,  it  does  not  provide  insight  into 
the  underlying  mechanisms  for  tensile  and/or  shear  deformations  in  filled  rubber.  This  is 
not  unexpected  since  model  (3.29)-(3.30)  is  based  on  pseudo-phenomenological  formulations. 
Hence,  a  different  approach  based  on  molecular  arguments  was  pursued  in  [4,  5,  6],  where 
the  ideas  of  these  models  are  based  on  those  of  Johnson  and  Stacer  in  [30] .  It  turns  out  that 
this  approach  leads  precisely  to  a  class  of  models  based  on  a  Boltzmann  formulation. 

A  polymer  material  undergoing  directional  deformation  is  modeled  in  [2]  in  which  polymer 
chains  are  treated  as  Rouse  type  strings  of  beads  interconnected  by  springs  (see  [45])  as 
depicted  in  the  left  plot  of  Figure  19.  The  model  permits  the  incorporation  of  many  impor¬ 
tant  physical  parameters  (such  as  temperature,  segment  bond  length,  internal  friction,  and 
segment  density)  in  the  overall  hysteretic  constitutive  relationship.  The  model  in  [2]  was 
based  on  the  assumption  that  the  materials  were  composed  of  two  virtual  compartments 
as  depicted  in  the  right  plot  of  Figure  19.  One  compartment  consists  of  a  constraining 


Piitiapperi  Portion  of  PC  molecule 


Figure  19:  (left):  Representation  of  vectors  for  a  bead-spring  polymer  molecule;  (right):  PC 
molecule  entrapped  by  the  surrounding  constraining  tube. 

tube  which  is  a  macroscopic  compartment  containing  both  CC  (chemically  cross-linked)  and 
PC  (physically  constrained)  molecules.  The  other  compartment  is  microscopic  in  nature  and 
consist  of  those  PC  molecules  aligned  with  the  direction  of  the  deformation.  These  molecules 
will  at  first  “stick”  to  the  constraining  tube  and  be  carried  along  with  its  motion,  but  will 
very  quickly  “slip”  and  begin  to  “relax”  back  to  a  configuration  of  lower  strain  energy.  In 
the  model  derivation  one  computes  the  contributions  of  both  “compartments”  to  the  over¬ 
all  stress  of  this  polymer  material  undergoing  deformations  to  obtain  the  constitutive  law, 
which  is  similar  to  that  developed  in  [4,  5]  and  has  the  general  form  of  Boltzmann  type 
model  (3.26),  even  though  the  kernel  is  not  of  convolution  type. 
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4  Example 


We  have  introduced  and  briefly  discussed  a  number  of  possible  constitutive  relationships 
that  represent  a  noncomprehensive  review  of  an  extensive  body  of  research  literature  on 
clastic/viscoelastic  materials.  We  conclude  by  presenting  an  example  to  illustrate  how  one 
might  apply  the  equations  of  motion  and  a  constitutive  relationship  to  model  a  particular 
situation.  In  our  computational  example,  we  will  examine  the  properties  of  a  one  dimensional 
column  of  dry  soil  and  study  the  movement  of  the  soil  in  response  to  a  sinusoidal  input 
at  the  surface.  At  first,  we  will  record  displacements  in  a  wave  moving  past  a  stationary 
observation  point  in  the  soil  column.  We  then  add  a  rigid  body  to  the  column  to  demonstrate 
the  changes  in  wave  propagation  one  might  expect  if  something  is  buried  in  the  soil.  We 
can  change  parameters,  such  as  soil  density,  and  understand  the  impact  these  changes  have 
on  the  displacements  at  the  observation  point.  In  all  cases,  we  focus  on  “seismic  P  waves” 
(longitudinal  waves)  propagating  downward  through  the  soil  away  from  the  source  of  the 
force  (an  impact)  located  at  the  ground  surface. 


4.1  Model  Description 

The  schematic  for  this  problem  is  given  on  the  left  side  of  Figure  20  where  the  wave  obser¬ 
vation  point  is  at  the  location  z  =  Zio  and  the  ground  surface  (as  well  as  the  source  of  the 
sinusoidal  input  force)  is  at  z  =  zp o.  The  right  side  of  the  figure  represents  the  idealized 
(1-dimensional)  mathematical  soil  column  with  input  source  and  observation  point  marked 
in  the  figure.  In  the  case  where  we  have  a  rigid  object  in  the  soil,  that  object  is  placed  at 
z  =  z±o  in  the  column  instead  of  an  observation  point. 
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Figure  20:  (left):  Schematic  of  problem;  (right):  1-dimensional  representation 


In  this  idealized  configuration,  we  assume  that  both  soil  and  target  are  uniform  in  x-  and 
^/-directions.  Based  on  the  discussion  in  Remark  3.4,  we  assume  that  soil  behaves  as  a 
Kelvin- Voigt  material  for  small  vibrations.  This  means  that  stress  components  in  soil  can 
be  expressed  as  sum  of  two  terms,  the  first  term  being  proportional  to  the  strain  (e)  and  the 
second  term  being  proportional  to  the  rate  of  change  (£)  of  strain. 
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Let  u(z,t )  denote  the  displacement  (in  units  m)  in  the  ^-direction  at  position  z  at  time  t. 
Then  in  the  situation  with  no  buried  object,  for  z  G  (zp o,  oo)  we  have 

Case  1:  Only  soil  present  in  the  column 


P(z) 


d2u(z, t ) 
d  t2 


d_ 

dz 


du(z,  t ) 
dz 


+  V(z) 


d2u(z, t ) \ 

dtdz  J  ’ 


(4.1) 


For  a  second  case,  we  assume  there  is  a  buried  rigid  object  with  its  center  of  mass  located  at 
z  =  z  io-  We  assume  that  the  target  is  homogeneous  in  the  z- direct  ion,  and  the  contacting 
surface  area  between  the  target  and  the  soil  under  the  target  is  the  same  as  the  one  between 
the  target  and  the  soil  above  the  target.  Since  the  target  is  a  rigid  body,  the  displacement 
of  upper  side  of  the  target  is  exactly  the  same  as  that  of  its  lower  side.  Hence,  in  the  1- 
dimensional  setting  we  can  treat  the  target  as  a  point  mass.  With  the  given  assumptions 
on  soil  and  target  we  can  visualize  the  soil-target-soil  as  two  thin  rods  connected  by  a  point 
mass  at  z  —  z io-  The  schematic  of  the  problem  is  illustrated  in  the  right  plot  of  Figure  20, 
where  in  this  case  there  is  a  point  mass  at  z  =  z io-  The  resulting  pair  of  equations  are 


Case  2:  A  rigid  target  present  in  the  soil  column 


p(z 


d2u(z, t ) 
dt2 


d  f  .  .  du(z,t )  .  .  d2u(z,t ) 


z  e  (ZpQ,Z10)  U  (zio,oo), 


M 


d2u(z10,t) 
dt 2 


=  S 


+  du(zt0,t)  +^d2u(zt0,t) 
MmoJ  -qZ  r  '7 Wo 


dtdz 


(4.2) 


-S 


_  du(z10,t)  _^d2u(z10lt) 

'^wio/  '  '/ Wio 


dtdz 


In  all  cases  p  denotes  the  density  (in  units  kg/m3)  of  soil,  k  is  the  clastic  modulus  (in 

units  — =  Pa)  of  soil  ,  and  r/  represents  the  damping  coefficient  (in  units  — — )  of  soil. 

For  the  second  equation  in  (4.2),  M  is  used  to  denote  the  mass  (in  units  kg)  of  the  target 
and  S  represents  the  surface  area  (in  units  m 2)  of  contact  between  the  target  and  the  soil 
under  (or  above)  the  target.  Though  the  model  will  treat  non-constant  and  piecewise-defined 
coefficients,  for  our  example  we  will  take  the  simple  case  where  the  parameters  are  constant 
values  in  the  soil  column.  In  general,  non-constant  and  piecewise-defined  coefficients  would 
allow  one  to  account  for  varying  physical  situations,  such  as  having  a  lower  density  above 
a  buried  object  than  below  the  buried  object.  For  initial  conditions  we  will  assume  zero 
displacement  and  zero  velocity  since  the  system  is  at  rest  initially.  These  conditions  are  then 
given  by 

du 

u(z,  0)  =  0,  twM)=0.  (4.3) 


dt 


The  boundary  condition  at  z  =  zp o  is  given  by 


.  .  du(z,t )  .  .d2u(z,t) 

k(z)  \  J  +  g{z N 


dz 


dtdz 


~zp0 


=  -fit), 


(4.4) 
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where  /  is  the  applied  external  force  in  units  N/m2.  On  the  surface,  the  normal  internal 
stress  is  balanced  with  the  input  force,  resulting  in  (4.4). 

In  order  to  numerically  solve  our  model  (4.1)  with  (4.3)  and  (4.4)  or  model  (4.2)  with 
(4.3)  and  (4.4),  it  is  convenient  to  have  a  finite  spatial  domain.  Since  we  only  care  about  the 
displacements  near  the  surface  (and  near  the  buried  object  in  the  second  case),  we  will  choose 
the  right  (lower)  boundary  £0o  to  be  sufficiently  far  away  from  the  target  so  that  no  energy 
will  reach  the  boundary  Zoo  during  the  time  frame  within  which  we  run  the  simulations.  This 
assumption  implies  that  we  can  set  up  any  type  of  boundary  condition  at  zoo-  For  simplicity, 
we  assume  that 

u(zoo,t)  =  0.  (4.5) 

Hence,  our  problem  (4.1)  with  (4.3)  and  (4.4),  and  problem  (4.2)- (4.4)  is  thus  defined  on  the 
finite  space  domain  [%,o,£oo]-  We  report  on  computations  for  the  model  using  a  standard 
finite  element  method. 


4.2  Simulation  Results 


We  report  the  results  of  some  of  our  simulations  with  the  equations  governing  this  one 
layer  problem,  observing  the  wave  form  through  time  at  location  z10  =  0.3048  m,  which 
is  approximately  one  foot  beneath  the  ground  surface  at  zp o  =  0.  The  value  for  the  far 
boundary  was  set  at  50  m,  and  no  reflections  from  the  far  boundary  were  observed  in  the 
calculations.  The  baseline  numerical  values  for  soil  density  and  elastic  modulus  are 

po  =  1800  kg/m3,  k0  =  2.04  x  108  Pa,  (4.6) 


where  the  values  were  estimated  from  [16].  For  the  value  of  damping  coefficient  of  soil  r/,  we 
will  assume  that  damping  is  frequency-  and  elastic  modulus-dependant,  via  the  formula 


2/3k 

w\/l  -I32 


(4.7) 


The  parameter  (3  is  called  the  damping  ratio,  which  is  related  to  the  energy  lost  between 
wave  peaks.  The  baseline  value  for  the  damping  ratio  was  set  at  /3 0  —  0.05,  which  again  was 
derived  using  results  from  [16].  For  the  carrier  frequency  in  our  simulations,  we  apply  the 
sinusoidal  input  function  at  a  frequency  of  to  =  4007T. 


We  depict  a  plot  of  the  sinusoidal  input  function  in  Figure  21.  Recall  that,  given  our 
coordinate  system,  positive  force  values  represent  downward  force  and  negative  force  values 
represent  upward  force.  Thus,  this  input  represents  impacting  the  ground  in  a  downward 
motion  and  then  the  ground  rebounding  with  equal  force.  We  can  thus  see  displacement 
return  to  the  baseline  in  all  of  the  following  figures  due  to  the  restoring  force  present  in  the 
input.  Since  the  displacement  returns  to  the  baseline,  we  can  also  infer  that  our  model  is 
linearly  dependent  on  the  input. 


42 


Forcing  Function 


Figure  21:  Sinusoidal  input  function. 


4.2.1  Results  for  Case  1:  Only  soil  present  in  the  column 

The  first  situations  we  examine  are  when  holding  one  parameter  (of  (k,  p))  constant  and 
increase  the  other  parameter.  The  results  are  depicted  in  the  panels  of  Figure  22.  In  the 
left  panel,  we  hold  density  p  constant  and  change  the  elastic  modulus  k.  As  the  modulus 
increases,  we  see  that  the  wave  begins  its  upward  slope  sooner,  meaning  the  wave  has  arrived 
at  the  observation  point  sooner  so  wave  speed  has  increased.  Also,  as  the  modulus  increases, 
we  see  less  displacement  overall.  In  the  case  where  density  increases,  we  see  that  we  also  get 
less  displacement  overall,  but  the  wave  speed  decreases.  We  clearly  see  the  model  demon¬ 
strates  that  soil  parameters  have  multiple  realistic  effects  on  overall  displacement  and  wave 
speed. 


x  -|  o-4  Displacement  around  z=0.3048m  (z=1  ft)  x  1 0-4  Displacement  around  z=0.3048m  (z=1  ft) 


Figure  22:  Wave  form  at  z\o  with  varying  soil  parameters  individually  (left:  variable  k, 
constant  p;  right:  constant  k,  variable  p ) 
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In  Figure  23,  we  examine  the  situation  where  we  increase  both  soil  parameters.  As  expected, 
the  effects  of  the  parameter  increases  combine  to  reduce  overall  displacement.  In  the  left 
pane  of  the  figure,  we  see  that  the  elastic  modulus  increases  more  rapidly  than  the  density 
and  the  wave  reaches  the  observation  point  sooner.  In  the  right  pane,  we  applied  the  same 
percentage  increases  to  the  parameters  and  so  the  wave  reaches  the  observation  point  at  the 
same  time  for  all  the  parameter  combinations.  Thus,  our  model  is  able  to  demonstrate  the 
effects  of  more  complex  parameter  changes  that  one  might  see  in  applications.  For  example 
one  might  input  multiple  impacts  at  the  same  site,  which  would  increase  both  soil  density 
and  stiffness  from  one  impact  to  the  next.  (This  in  fact  was  realized  in  field  experiments  by 
scientific  colleagues  when  testing  for  repeatability  of  responses  to  interrogation  impacts.) 


Figure  23:  Wave  form  at  Z\q  with  both  soil  parameters  k  and  p  varying. 


4.2.2  Results  for  Case  2:  Rigid  body  present  in  the  soil  column 

In  this  section  we  examine  the  case  where  we  have  included  a  rigid  body  located  at  zio, 
modeled  as  a  point  mass  at  the  same  point  as  our  previous  observation  point.  In  the  figures, 
this  location  is  represented  by  a  vertical  dashed  line.  For  all  the  simulation  results  presented 
in  this  section,  the  values  for  the  soil  density,  elastic  modulus  and  damping  ratio  are  chosen 
to  be  the  baseline  values  specified  earlier,  with  the  value  for  the  mass  of  the  target  given  by 
M  =  34.2671  kg ,  and  the  contacting  surface  area  taken  as  S  —  0.1580  nr? . 

In  the  upper  left  pane  of  Figure  24,  we  can  see  that  the  sinusoidal  force  has  begun  imparting 
displacement  in  the  soil  but  this  displacement  has  not  yet  reached  the  target.  In  the  upper 
right  pane,  the  displacement  has  impacted  the  target.  In  this  simulation,  the  target  imparts 
much  of  the  energy  to  the  soil  below  it.  In  the  lower  left  pane  of  Figure  24,  we  see  that  most 
of  the  energy  has  passed  through  the  target  and  is  deeper  into  the  soil.  However,  looking 
at  the  domain  between  z  =  0  and  the  dashed  target  line,  we  see  the  remnant  energy  that 
was  reflected  by  the  target  back  toward  the  surface.  In  the  lower  right  pane  in  the  figure, 
we  see  that  this  energy  has  bounced  off  the  soil  and  impacted  the  target  again,  once  more 
passing  some  energy  through  the  target  (as  seen  in  the  small  amount  of  extra  displacement 
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Wave  form  in  z-domain,  at  time  t=0.00067831  x  ^  g-4  Wave  form  in  z-domain,  at  time  t=0.0020486 


Wave  form  in  z-domain,  at  time  t=0.0054744  x  -|  g-4  Wave  form  in  z-domain,  at  time  t=0.0068447 


Figure  24:  Wave  form  in  the  soil  column  at  various  times  (the  buried  object  location  is 
represented  by  vertical  dashed  line). 


to  the  right  of  the  dashed  line)  and  having  some  energy  reflected.  The  basic  one  dimensional 
model  we  have  used  can  effectively  model  the  presence  of  a  rigid  body  in  a  column  of  soil, 
in  particular  the  behavior  of  partial  reflecting  and  transmitting  of  energy  when  the  object  is 
impacted. 

This  numerical  example  has  demonstrated  the  ability  of  a  one  dimensional  model  to  capture 
some  of  the  salient  features  of  wave  propagation  in  the  soil  medium,  including  sensitivity  to 
soil  structure  and  to  the  presence  of  rigid  objects  in  the  soil.  Additional  features  that  would 
require  a  two  or  three  dimensional  model  might  be  modeling  the  presence  of  more  than  one 
type  of  body  wave  (e.g.,  shear  and  compressional)  as  well  as  modeling  surface  waves.  In  such 
a  higher  dimensional  setting  we  could  also  take  into  account  more  complicated  buried  object 
geometries.  Ultimately,  the  one  dimensional  model  still  captures  much  of  the  basic  dynamics 
of  elasticity  in  soil  and  is  therefore  useful  in  predicting  outcomes  to  physical  experiments 
with  some  degree  of  fidelity. 
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